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Context. Abridged. Many stars are members of binary systems. During early phases when the stars are surrounded by a discs, the binary orbit 
and disc midplane may be mutually inclined. The discs around T Tauri stars will become mildly warped and undergo solid body precession 
around the angular momentum vector of the binary system. It is unclear how planetesimals embedded in such a disc will evolve and affect 
planet formation. 

Aims. We investigate the dynamics of planetesimals embedded in gaseous protoplanetary discs that are perturbed by a binary companion on a 
circular, inclined orbit. We examine the collisional velocities of the planetesimals to determine when planetesimal growth through accretion can 
occur, instead of disruption. We vary the binary inclination, y F , binary separation, D, disc mass, M d , and planetesimal radius s t . Our standard 
model has D = 60 AU, jf = 45°, and a disc mass equivalent to that of the minimum mass solar nebula model. 

Methods. We use a 3D hydrodynamics code to model the disc. Planetesimals are test particles which experience gas drag, the gravitational 
force of the disc, the companion star gravity. Planetesimal orbital crossing events are detected and used to estimate collisional velocities. 
Results. For binary systems with modest inclination (y F = 25°), disc gravity prevents planetesimal orbits from undergoing strong differential 
nodal precession (which they would do in the absence of the disc), and forces the planetesimals to precess with the disc on average. For 
planetesimals of different size the orbit planes become modestly mutually inclined, leading to collisional velocities that inhibit planetesimal 
growth. For larger inclinations (y F = 45°), the Kozai effect operates, leading to destructively large relative velocities. 

Conclusions. We conclude that planet formation via planetesimal accretion is difficult in an inclined binary system with parameters similar to 
those considered in this paper. For highly inclined systems in which the Kozai effect switches on, the prospects for forming planets are very 
remote indeed. 
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1. Introduction 

Of the extrasolar planets detected so far, over 20 percent are 
found to orbit one component of a multiple/binary st ar system 
( Des idera & Barbi eril 2007 : EggenbergeretaL 2007 1. Planet 
formation in binary systems can represent a particular chal- 
lenge, as each stage of the formation process can be affected 
by the gravitational perturbation of the binary companion. One 
crucial stage that is particularly sensitive to the presence of the 
companion star is the accumulation stage of kilometre-sized 
planetesimals. 

The fundamental parameter that controls the efficiency of 
planetary growth by accretion of planetesimals is the rela- 
tive velocity between impacting objects. This must be lower 
than the escape velocity from the accreti ng objects in order 



for e fficient runaway accretion to occur ( Wet herill & Stewart 
If the external gravitational perturbation by the bi- 



nary companion excites relative velocities that exceed the es- 
cape velocity, runway accretion is prevented and growth re- 



mains slow. Furthermore, if the relative velocity is excited 
beyond the threshold velocity at which erosion dominates 
accretion, planetesimal growth potentially ceases altogether 
feenz& Asphaugfll999l) . 



The possible effect of a stellar companion on the perturbed 
velocity distribution of planetesimals has been explored in sev- 
eral previous publications, where most have considered con- 
figurations in which the binary orbit is eccentric and coplanar 
with t h e disc midplane ( Hepp enheime rl 1 978t Ma rzari & Scholl 
| 2000[ IThebault. Marzari & Scholll 120061: IPaardekooper et all 
2008). The main conclusion of these studies was that the cou- 
pled effects of secular perturbations of the binary companion, 
and friction due to gas drag arising from the protoplanetary gas 
disc, induce forced eccentricities and a size-dependent phasing 
of pericentres. This leads to relatively modest collision veloc- 
ities dominated by the keplerian shear for same-sized bodies, 
but high relative velocities for different sized planetesimals. As 
a consequence, binary systems with separations ~ 50 AU may 
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have a strong inhibiting effect on accretion within a swarm of 
colliding km-size planetesimals. 

These studies neglected the effect of disc gravity acting on 
the planetesimals, which can affect the details of the results 
when the disc becomes eccentric through interaction with a bi- 
nary system on an eccentric orbit. But, the main conclusion that 
planetesimals of different sizes experience large collisional ve- 
lociti es due to differential p hasing of their pericentres remains 
valid dKlev & Nelsonll2008l) . 

These previous studies assumed that the orbital plane of 
the binary co mpan i on is coplanar with the disc midplane. 
According to iHald (119941) the assumption of coplanarity or 
modest inclination may be reasonable for binary separations 
below ~ 40 AU, but systems with larger separations appear to 
have their orbital planes randomly distributed. The examination 
of planetesimal dynamics in non- coplanar configurations has 
received attention only recently dMarzari. Thebault & Scholll 
2009:Xie &~Zhoul2 009). The former authors show that, due to 
the semi-major axis dependence of the nodal precession rate, 
the nodal lines of the planetesimals become progressively ran- 
domized. This may lead to the dispersion of the planetesimal 
disk that expands into a cloud of bodies surrounding the star. 
The latter authors considered the effect of gas drag in an sys- 
tem with a modestly inclined binary, and showed that in addi- 
tion to inducing size-dependent phasing of pericentres, gas drag 
also introduces phasing of nodal lines. This leads to a situation 
is which pl anetesimals of diffe rent size occupy different or- 
bital planes. IXie & Zhoul (120091) suggest that this may provide 
a favourable channel for planetesimal growth as low velocity 
collisions between similar sized objects become more frequent 
than high velocity collisions between different sized objects. 
It should be noted, however, that these authors neglected the 
effect of disc gravity on the planetesimals. Although the inclu- 
sion of disc gravity modifies the dynamics of planetesimals in 
fully coplanar systems, the changes are not dramatic. In non- 
coplanar systems, however, where the disc and planetesimals 
may precess at different rates around the binary angular mo- 
mentum vector, the deviation of the planetesimals from the disc 
plane results in the disc gravity providing a strong restoring 
force that can modify the dynamics in an important qualitative 
sense. 

Th e r ecent simula t ions by iMarzari. Thebault & Scholl 
(2009) and Xie & Zhou (2009) either ignored the gaseous pro- 
toplanetary disc, or treated it as a static object which did not 
evolve in the gravitational field of the inclined companion. It 
is known, however, that a gaseous disc orbiting under the in- 
fluence of a close binary companion on an inclined orbit will 
develop a mild warp and precess as a solid body around the or- 
bital angular momentum vector of the binary system if bending 
waves can propagate efficiently (Paoaloizou & Prinde 1983; 



Panaloizou & Lin 1995 
1995; Lubow& Ogilvie 



Odlvie 2000; Papaloizou & Terquem 



20001) . T his process has been s tudied 



using numerical simulations by lLarwood et al.l (119961) . who 



used Smo oth Particle Hydrodynam ic calculations, and more re- 
cently by Fragner & Nelson ( 2010i) . who performed 3D simu- 
lations using a grid-based hydrodynamics code . The condition 
for bending waves to propaga te is that a < H/R, where a is 
the usual viscosity parameter (IShakura & Sunvaevll 19731) and 



H/R is the disc aspect ratio. In protostellar discs, it is estimated 
H/R ^ 0.05 and a ^ 10~ 3 - 10~ 2 , so we expect such discs to 
evolve similarly to the models presented in this paper. In addi- 
tion to generating a mild warp, and causing the disc to precess, 
the binary can also tidally truncate the disc at its outer edge, 
where the tidal truncation radius of the disc for a binary system 
of unit mas s ratio is typically ~ D/3, where D is the binary 



iry 

separation dArtvmowicz & Lubow 1994 : lLarwood et aD ll996: 
Fragner & NelsonlbOlOh . 



An observational e xample of a disc in a misaligned bi- 
nary system is HK Tau ( Stanelfe ldt et alJ ll998). where the bi- 
nary system and disc have been imaged directly. More indirect 
evidence for precessing discs in close binary systems comes 
from observations of precessing jets in star forming regions 
dTerquemetal.lll999h . 



In this paper we investigate the dynamics of planetesimals 
embedded in protoplanetary disc models in inclined binary sys- 
tems with circular orbits. Although most binaries are on eccen- 
tric orbits, we chose the case of a circular orbit to allow us to 
focus on the effects inclination. In contrast to previous work, 
we solve for the disc evolution in the gravitational field of the 
inclined binary companion, and include the effects of the disc 
gravity acting on the planetesimals. We consider two different 
values of the inclination between the binary orbit plane and disc 
midplane, jp - 25° and jp = 45°, planetesimals of size 100 m, 
1 km and 10 km, and binary separations between 60 and 120 
AU. The disc outer radius is 18 AU, consistent with the tidal 
truncation radius for a binary separation D = 60 AU. 

Due to the complexity of our model, and the associated 
computational expense of running the simulations, we are un- 
able to consider large numbers of planetesimals, and so we 
are forced to use an approximate method for determining their 
collisional velocities based on examining the moments when 
the osculating orbits of the planetesimals intersect. Using this 
approximation, we find that in general the binary companion 
causes large planetesimal collision velocities to be generated, 
largely because the orbit planes of the planetesimals become 
mutually inclined. For the larger value of the inclination, we 
find that the Kozai mechanism can switch on, leading to the 
generation of large orbital eccentricities for the planetesimals, 
and therefore very large collision velocities. These results in- 
dicate that planet formation via the accumulation of planetesi- 
mals will be difficult in binary systems with parameters similar 
to those we consider in this paper. 

This paper is organised as follows. In Sect. 2 we present 
the basic equations and in Sect. 3 we discuss the numerical 
techniques. In Sect. 4 we investigate numerically the effect of 
disc gravity on the evolution of inclined planetesimal orbits in 
the absence of a binary companion. In Sect. 5 we examine the 
dynamics of planetesimals that are embedded in a disc which 
is inclined initially with respect to the binary plane, and ex- 
amine the effect of varying the inclination angle jp. In section 
Sect. 6 we present calculations for different disc masses and bi- 
nary separations and examine under which condition the Kozai 
effect can be suppressed. Finally, we discuss our results and 
draw conclusions in Sect. 7. 
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2. Basic equations 

The equations of hydrodynamic s for a viscous disc th at we 
solve in this paper are given in iFragner & Nelsonl fcoioh . The 
disc evolves under the effects of pressure, viscosity and the 
gravitational forces due to the binary companion and central 
star. We employ a large number of variables in this paper to 
describe the results, and we tabulate these in table Q]for ease of 
reference. 

We solve the equations in a frame that precesses around 
the angular momentum vector of the binary orbit, since it is 
expected that the disc models we consider in this paper will 
undergo uniform precession due to interaction with the binary 
companion. In this frame the disc midplane stays close to the 
equatorial plane of the spherical polar grid that we adopted in 
the simulations, provided that the precession frequency of the 



pr 

frame, Q F , is chosen according to (Fra gner & Nelson]|2010l) : 



£1 F 3R 3 M B 



(i) 



where Cld(R) is the orbital angular velocity of the disc at its 
outer radius R, D is the binary separation, and Mb are the 
masses of the primary and secondary star, respectively, and jf 
is the relative inclination between the disc midplane and the 
binary orbital plane. 

The binary companion is held on a fixed circular orbit with 
separatio n D. Its position vector, D , in the precessing frame is 
given by ( Fragner & Nelsonll2010l) : 



D = D 



COS ([(l>b - 0.f]t) 

sin ([cj b - Qf]f) cos(y F ) 
sin([w B - Q F ]f))sin (y F ) ) 



(2) 



where o>b is the binary angular frequency measured in the non- 
precessing binary frame. Thus an observer moving with the 
precessing frame sees an increased binary frequency a>B - Of 
due to the retrograde precession of the frame (i.e. Q F is nega- 
tive). 

The planetesimals we model do not mutually interact. They 
feel the gravitational force of the primary star, secondary star 
and disc, the drag force due to the disc, Coriolis and centrifugal 
forces due to the precession of the frame and indirect forces that 
arise because we centre our coordinate system on the primary 
star. The equations describing the evolution of the planetesi- 
mals are therefore: 



d 2 r, 

dt 1 



GM+ 



•i-GM B p °3 -G f dm(r) 
|r» - DP J disc 



fo - r) 
h-rp 



— - 2£l F x Vj - £If x (£l F x rO 



GM B 



D-G f i 

J disc 



r 

dm(r) — 



(3) 



where G is the gravitational constant, r,, V; and M, are the po- 
sition and velocity vector and mass of planetesimal i. The first 
and second terms are the gravitational acceleration due to the 
central and companion star, respectively, and the third repre- 
sents the gravitational acceleration due to the disc. The fifth 
and sixth terms represent Coriolis and centrifugal forces that 



arise because the coordinate system pr ecesses around a vector 
£If , given by (Fragner & NelsorfeOlOh : 



£l F = Q.p 





- sin(y F ) 
cos(y F ) ) 



(4) 



The last two terms are indirect terms that account for the accel- 
eration of the central star by the companion and disc, respec- 
tively . The gas drag force Fo can be written (Weidenschilling 
1977h: 



2 1 

F D = -ConSiP-Wi - v|(Vj - v) 



(5) 



Here s, is the radius of planetesimal i, p and v are the gas disc 
density and velocity. The value of the drag coefficient Co de- 
pends on the Reynolds number 



K. 



2si\Vi ■ 



rmol 



(6) 



where the molecular viscosity, v mo ; = Amc s , Am is the mean 
free path of the gas molecules, and c s is the sound speed. Co 
takes values: 



24. -r; 1 %, < 1 
C D = 24. ft; 6 1 < R e < 800 
0.44 % e > 800 

3. Numerical method 



(7) 



The hydrodynamic disc equations are inte grated using the 



grid-b ased hydrodynamics code NIRVANA (IZiegler & Yorke 



1997), adapted to solve the equations in a precessing reference 
frame. This code uses operator-splitting, and the advection rou- 
tine u ses a second-or der accurate mono tonic transport algo- 
rithm (Ivan Leerlll977l) . The planetesimal orbits are integrated 
using the leap-frog integrator. 

3.1. Units 

The equations are integrated in dimensionless form, where we 
choose our unit of length to be the radius of the disc inner edge. 
The unit of mass is that of the central star. The unit of time used 
in the code corresponds to f2^(l) _1 (where 0#(1) is the keple- 
rian angular frequency at radius R — 1), and the gravitational 
constant G = 1. When discussing simulation results we will 
refer to a time unit that corresponds to Pj = 27t/Qa:(10), which 
is approximately one orbit at the outer edge of the disc. In the 
sections of this paper which describe the results of simulations 
we refer to evolution times in units of "orbits", where an orbit 
should be understood as being a time interval equal to Pj. We 
scale our unit of length to 2 AU in physical units, and assume 
that the central star has one solar mass, so that one orbital pe- 
riod at the outer disc edge corresponds to a time of Pa = 89.44 
years. Inclination and precession angles are displayed in units 
of radians. 
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3.2. Initial and boundary conditions 

The disc mode l we use is model 1 described in 
Fragner & Nelson (1201 rl . and interested readers are re- 
ferred to that paper for a general description of disc evolution 
in misaligned binary systems. The disc model extends from 
1 - 9 in code units, corresponding to 2 - 18 AU in physical 
units, has a height to radius ratio h - 0,05, and a dimensionles s 
viscosity parameter a = 0.025 dShakura & Sunvaevlll973l) . 
During its evolution, the disc model develops a very modest 
warp (the variation in inclination across its radial extent is 
less than one degree), and it precesses as a rigid body around 
the angular momentum vector of the binary system with 
a frequency approximately equal to Sip given by Eq. (Q~|i. 
The disc exhibits spiral density waves which are excited by 
the companion, and does not appear to become noticeably 
eccentric during its evolution - unlike discs which are per- 
turbed by lower mass binary companions (Kiev et al. 2008J), or 
companions on eccentric orbits (IKlev & Nel son 2008). 

In our standard model we normalise the disc mass so that it 
would contain 0.015 M if extended out to a radius of 40 AU 
(this being very similar to the mass co ntained in the minimum 
mass solar nebula model dHavashilll98lR although the model 
we use nominally only goes out to ~ 18 AU. The actual disc 
mass is M d = M = 4.52 x 10~ 3 M G , where M denotes our 
standard disc mass. In physical units the gas disc density is 
given by: 



p(r,0) = 7.1 



sin 



0p- 



g cm 



-3 



where is the meridional angle measured relative to the nor- 
mal to the disc midplane. The binary companion is held on 
a circular inclined orbit with constant separation D — 30 
(= 60 AU). We note that mos t binaries are on eccentric or- 
bits (iDuquennov & Mavorlll99ll) . but we restrict ourselves to 
circular inclined orbits in this study so as to isolate effects re- 
lating to the inclination. Note that we neglect the disc gravity 
acting on the binary companion, so the binary orbit is fixed. At 
the beginning of each simulation, the companion mass is in- 
creased linearly until it reaches its final mass Mb = 1M* after 
4 orbits. Its angular frequency u>b is increased accordingly to 
be consistent with a stationary orbit. 

Periodic boundary conditions were applied in the azimuthal 
direction. At all other boundaries standard stress-free, outflow 
conditions were employed. 

3.3. Planetesimal set-up 

The mass of the planetesimals is given by M, = ins^ps, and 
we choose a solid density of p s = 2 g cm' 3 . They are initially 
embedded within the unperturbed disc model on circular, kep- 
lerian orbits which are coplanar with the disc midplane. As the 
disc does not become noticeably warped or eccentric during its 
evolution, it is not necessary to pre-evolve the disc so that it 
achieves a steady state structure prior to inserting the planetes- 
imals. 

In the following discussion we will characterise the plan- 
etesimal evolution using their orbital elements in the fixed bi- 



nary frame, and these are calculated from the positions and ve- 
locities that the code outputs in the precessing frame. Hence, 
we first have to transform the velocities and positions from the 
precessing frame into the binary frame, in which the angular 
momentum vector of the binary Jg is parallel to the unit vector 
e 3 . 

As we consider two reference frames in this paper, one 
which precesses with the disc, and the fixed frame based on the 
binary system, we denote vectors and coordinate values in the 
precessing frame using the hat-symbol (e.g. f, 9). Vectors and 
coordinate values in the fixed frame are denoted without the 
hat-symbol (e.g. r, 6). The transformation of the position and 
velocity vectors f , and v, defined in the precessing frame into 
vectors r, and v, defined in the no n-precessing binary frame is 
given by ( Fragner & Nelsonll2010l) : 



r, = R x \y F )Rz(^Ft)h 
d 



at 



(8) 



where Rz and Rx are rotation matrices around the z and x axes, 
respectively. The last term accounts for precessional velocities. 
Note that in a strict sense, we should replace Q.pt by J Qpdt 
since the precession frequency is increased during the first 4 
orbits of the simulations. For simplicity, however, we keep this 
notation and understand that this replacement should be made. 
The velocity and position vectors can be used to calculate or- 
bital elements in the binary frame. These are denoted a,, e,, Q,, 
a,-, o>i, fi for the semi-major axis, eccentricity, longitude of as- 
cending node, inclination with respect to the binary plane, lon- 
gitude of pericentre and true anomaly of particle i, respectively. 
Since the planetesimals are initialised to be coplanar with the 
disc midplane, and the equatorial plane of the spherical com- 
putation grid, it follows that a, = jp and £2; = n at t = 0. 

Additionally, we are interested in the inclination of the 
planetesimals with respect to the local disc midplane, which 
we will denote by the symbol 5: 



cos (Si) 



|L(n)l Ud 

sin (a,) sin (av/) cos (£2, - Q.^) + cos (ay) cos (af) (9) 



where j, is the specific angular momentum vector of particle 
i, and L(r,) is the angular momentum vector of the disc an- 
nulus which has the same distance from the central star, r,-, as 
planetesimal i at the current time. In this way, we measure the 
inclination of the planetesimal with respect to the local disc, 
which we will describe as the relative inclination from now on. 
The symbols a d and denote the inclination and nodal pre- 
cession angles, respectively, of the local gas disc annulus with 
respect to the binary orbital plane. Because the disc is rigidly 
precessing (and not strongly warped), and because of our ac- 
curate choice of Q.f, the disc midplane stays very close to the 
equatorial plane of the computational grid, and we have ap- 
proximately that = jp and £2^ — n + Clpt during the sim- 
ulations. If planetesimals stay close to the disc midplane (such 
that |Q, - £2,/| <s 1, Iff; - ad\ ^ 1), Eq. (|9]l can be approximated 
by: 



8] - (Q ; - a d ) 2 sin 2 («,■) + (a,- - a d y. 



(10) 
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Note that such an approximation is generally valid for two or- 
bits, whose orbital parameters (£2,, ai) are very similar. 

It will be instructive to understand which contributions 
(gas drag, disc gravity, binary companion) are responsible for 
changing the orbital elements of the planetesimals. In order 
to measure this numerically, each force contribution is trans- 
formed into the binary frame, where we take the sum of the di- 
rect and indirect parts when considering the gravity of the disc 
or companion. We then calculate the radial Fr, tangential Fj 
and normal F^ components with respect to the planetesimal 
orbit in the binary frame for each of the accelerations. These 
can be used to calculate the cha nges of osculating orbital ele- 
ments dMurrav & Dermottl 19991) . Here we will only state some 
of them that will be used later for the discussion of results (not- 
ing that a dotted quantity denotes its time derivative): 



di = Fn ■ — cos + fi) 
ji 



Q - F N ■ — 



Iji - (e 3 .j/)e 3 | 



sin (ait + f). 



(ID 
(12) 



For the change of relative inclination, differentiating Eq. (0 
with respect to time gives: 

- sin (di) ■ 6i = cos (ai) sin (ad) [an cos (£2,- - £2</) - a D ] 
+ sin (a,-) cos (ad) [«d cos (£2,- - Qj) - a,] 
+ sin (a,-) sin (ad) sin (£2, - Cld)(Cld - £2,). (13) 

When calculating the change of relative inclination Si we note 
that changes in the disc inclination or nodal precession can only 
be caused by the binary companion. Hence when considering 
the change due to drag force or disc gravity we set ad = and 
Q rf = in Eq. (13). 

3.4. Collision detection 

In order to obtain collisional velocity distributions that are sta- 
tistically relevant, it is necessary to accumulate data on a large 
number of direct collisions between planetesimals. Because we 
include the effect of disc gravity acting on the planetesimals, 
which incurs a large computational overhead, we are only able 
to integrate 50 planetesimals for each size that we consider. 
This means that we need to use an approximate method for es- 
timating the collision velocities between planetesimals, instead 
of detecting direct collisions between them. 

The approximation that we adopt in this work is to treat the 
planetesimals not as individual particles, but rather as repre- 
sentatives of their orbits, which we conceive of as being highly 
populated. These orbits are assumed to have a circular cross 
section with a radius Ar = 2 • 10~ 4 AU in physical units (note 
that this value has been used for the inflated radii of planetes- 
imals in Kie & Zhoul (2009), where a direct collision detec- 
tion method was used). In order to estimate impact velocities 
for colliding planetesimals, we detect the moments when two 
orbits defined by their osculating elements intersect, and es- 
timate the velocity of impact that planetesimals at the point 
of intersection would have. When running simulations which 
consider the general dynamics of planetesimals in misaligned 
binary systems, we distribute the planetesimals with a large 



range of semi-major axes within the protoplanetary disc ini- 
tially. When running simulations which are designed to cal- 
culate the collisional velocities, however, we perform separate 
runs in which the initial planetesimal orbits are distributed with 
a narrow range of semi-major axes (Aa = 10~ 3 AU). This is to 
maximise the number of orbit crossings, and hence improve our 
collision statistics. 

The condition for two orbits represented by particles i and 
j to cross is given by the vector equation: 



r,(0/, Cli, ai, at, e t , u)j) = r j(<p j, Clj, a h a h e h ojj) 
where r,- is given by: 



(14) 



r;(0i, a u e u a>i) 



' cos (£2,) cos ((pi) - sin (£2,) sin ((pi) cos (ai) " 
sin (£2j) cos ((pi) + cos (£2,) sin ((pi) cos (ai) 
sin ((pi) sin (ai) 



(15) 



with 



n(<f>i, a u e u u>i) = 



a t (\ -eh 



1 + e,- cos ((pi - LOj) 1 + e ; cos ( fi) ' 



(16) 



and similarly for the orbit represented by particle j. These vec- 
tor equations involve orbital elements which are defined in the 
fixed binary frame. The angles (pi and <pj in Eq.(fT4li are the an- 
gular distances of the orbit crossing point to the nodal lines of 
orbits ; and j, respectively, where the nodal lines are located 
in the binary orbit plane. For general crossing orbits, where 
the semi-major axes and eccentricities of the two orbits differ, 
it is not possible to solve Eq. ( fT4b directly. For circular orbits 
with the same semi-major axis (r,- = r X however, we can solve 
Eq. ( TBI for (pi and (pj, giving the values of these angles at the 
two points of intersection for orbits i and j: ((pj\, <pn) and (<pfi, 
(pi2). These angle are defined by the expressions 



tan (0/0 



sin (£2, - £2j) 



(17) 



cos (af) cos (£2,- - £2^) - sin (af)l tan (ai) 
0/2 = <pj\ + n 
cos ((pii ) = cos (<pj] ) cos (£2/ - Q.j) 

+ sin ((pji) cos (aj) sin (£2,- - £2y) 

sin(0 7 i)sin(a ; ) 

sin(0;i) = —. 

sin (ai) 

The solution for <pa is obtained by using (pj2 instead of (pj\ in 
the last two equations. 

The above solutions have been derived under the assump- 
tion that the orbits are circular with the same semi-major axis. 
If they have a finite eccentricity and different semi-major axes, 
we still assume that the point of closest approach is at these 
longitudes^ ((pn, <pj\) and ((pa, <pj2). Hence we define orbital 
crossing or collision, if the following condition is fulfilled: 

\ri((pn,ai,ei,(i)i)- rj((p n ,aj,ej,Uj)\ < Ar (18) 

where ri(<pn,ai,ei,ui) is given by Eq. ( fTSI i. and this condition 
is also checked for the other potential crossing point ((pa, <pji)- 



1 Note that this simplification was used to speed up the collision test 
in the simulations. Full numerical evaluation of the closest approach 
point of the two orbits gave good agreement with this assumption. 



6 



M.M. Fragner, R.P. Nelson & W. Kley: On collisional growth of planetesimals in misaligned binary systems 



Variable 


Variable description 


Jf 


Inclination between disc midplane and binary orbit plane 


D 


Binary separation 


(o B 


Binary orbit frequency 


fl F 


Precession frequency of the precessing frame and disc 


Pd 


Orbital period at 20 AU 


Mi 


The disc mass used in a simulation 


M 


Our nominal disc mass (M = 4.52 x 1(T 3 M G ) 


h = H/R 


Disc aspect ratio 


a 


Viscosity parameter 


P 


Gas density 


Si 


Planetesimal physical radius 


X 


Vectors and coordinate quantities in the precessing frame denoted with the 'hat' symbol 


X 


Vectors and coordinate quantities in the fixed binary frame denoted without the 'hat' symbol 


Oj 


Semi-major axis of planetesimal i 


e; 


Eccentricity of planetesimal i 


ft 


True anomaly of planetesimal i 


or; 


Inclination angle between binary orbit plane and orbit plane of planetesimal (' 


OJj 


Longitude of pericentre of planetesimal i (also referred to as the apsidal precession angle) 


fi, 


Longitude of ascending node of planetesimal i (also referred to as the nodal precession angle) 


Si 


Inclination of planetesimal i relative to the local disc 


a d 


Inclination of local disc annulus relative to the binary orbit plane 




Precession angle of local disc 


at 


Inclination angle between binary orbit plane and orbit plane of planetesimal i 



Table 1. Table of variable names. 



After every one hundred time steps the condition ([18} is 
checked for the 1 1 175 particle pairs that arise from the 50 par- 
ticles integrated for each size. Since this gives a total on the or- 
der of 10 7 pairs for the simulated time intervals we considered, 
we obtain quite a large number of pairs for which condition 
(fj~8b is fulfilled. Hence we are able to obtain good statistics on 
collisional velocities despite the relatively low number of plan- 
etesimals modelled. 

For orbit i the velocity at the first orbit crossing location is: 

d 

V;i = — r,(0a , n„ a,-, <?,-, «,) = v„r,- + v^<f>i 



\|a,(l-ef) L J 

(19) 

and likewise for orbit j. Here f , is the particle unit vector calcu- 
lated from Eq. ([TBl l and 0,- is the unit vector in the (f> direction. 
Likewise v„ and v^,- are the velocity components in the r and (p 
directions. Note that this expression does not account for pre- 
cessional velocities. However, since they do not contribute to 
the relative velocities at orbital crossing points (where r,- = rj) 
they can be safely omitted. The relative velocity at the first 
crossing location is then simply Av = |v,i — V/i |, and likewise for 
the second location (ftp). Later we will present averages of 
collisional velocities calculated in this way. 

3.4.1 . Validity of the orbit crossing technique 

An important question to address is under what conditions does 
the orbit crossing technique that we have described above agree 



with the direct detection of collisions when calculating average 
collision velocities within a planetesimal swarm. 

The first point to note is that the method we use to calcu- 
late the point of closest approach between two orbits assumes 
that this occurs along the line of intersection between the two 
orbit planes. For orbits which are mutually inclined with re- 
spect to one another, this assumption is valid. But it obviously 
breaks down for coplanar orbit planes, and our method reports 
the wrong points of intersection in this case. So our method of 
detecting orbit crossing is only strictly valid for mutually in- 
clined orbits. 

In order to address the more general question of whether 
the orbit crossing method gives collisional velocities which are 
similar to the direct collision detection method, we have run a 
number of pure N-body simulations. The first test we ran used 
a narrow ring of particles centred around 10 AU with static 
orbital elements (no binary companion was included) that were 
very similar to those reported in Fig. 7 of this paper at 60 orbits. 
The direct detection simulation used 1000 particles and was run 
for approximately 2000 local orbits. The orbit crossing simula- 
tion run used 50 particles and was also run for approximately 
2000 local orbits. We found that the mean collision velocities 
reported by each method agreed to within approximately 30 
%. The second test was similar, but included a binary compan- 
ion inclined by 25° to the ring of planetesimals, which were 
initially placed on circular orbits. This was run for 2000 local 
orbits and again good agreement was found between the two 
collision detection methods, with the mean collision velocity 
reported by the orbit crossing method being approximately 40 
% larger than that reported by the direct detection method. We 
conclude that the orbit crossing method gives fairly reliable re- 
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suits for the mean collision velocity in a planetesimal swarm 
provided that the dominant contribution to the relative veloci- 
ties arises because of differential nodal precession. 

3.4.2. Analytical estimate of collisional velocities 

To gain a better understanding of the contributions that control 
the collisional velocities, we now derive an analytical estimate. 
Consider two particle orbits A and B, which have very similar 
orbital elements (i.e. «a = o-b, e B — <?a + Ae, Q^ = Qa + AQ, 
as — c?a + Aor, <x>b - <jl>a + Aw), where the A-quantities are 
assumed to be very small. The above orbital elements are mea- 
sured with respect to the binary orbit plane. Consider now a 
coordinate system that is coplanar with the orbit of particle 
A. The velocity of particle A is then given by Eq. (TT9b - i.e. 
va = v r A?A +V0a$a • With respect to the x—y plane of this coordi- 
nate system, the orbit of particle B will be inclined by an angle 
Sab, because its plane is slightly different due to the quantities 
Aa and AQ. But since these differences are very small by as- 
sumption, the relative inclination Sab between the two particle 
orbits is small. Note also that in such a coordinate system we 
can always choose Qa such that Q^ = 0. Hence from Eq. (15[ 
it follows that: 



r B = ?a + dr with dr = 






sin (4> b )Sab 



(20) 



where Eq. (l20l has been derived using the assumption that the 
mutual inclination Sab ^ 1. Hence the orbit crossing point 
(ta = Tb) corresponds to 4>b = 0, from which it follows that: 



dr = and d(p = 



( \ 


Sab , 



(21) 



The velocity vector of orbit B at the orbit crossing point is thus: 



v B = (v r A + dv r )f A + (V^A + ^V )(0a + d</>) 



(22) 



To first order we neglect the term dv^dtp, and the relative ve- 
locity between the two orbits at orbital crossing is: 

Av = \ B - Va = dv r f A + dv^A + v Ad0. (23) 

pA-d(p = 0, the relative velocity 



Because Ta-4>a = ?A-d0 = 
becomes: 



Av 2 = dv 2 r + dv 2 + v 2 S\ B = dv 2 r + dv 2 + v\s\ B . 



(24) 



In the second equality we can replace v^a by the keplerian ve- 
locity vk to this order. From Eq. ([19} the radial and azimuthal 
velocity differences are to first order: 

dv r = v r (e + Ae, u + Aw) - v r (e, a>) 



GM+ 



\a(l -e 2 ) 
dv$ = v^(e + Ae, u> + Aa>) - v<p(e, of) 



(Ae sin ((/> — cS) — eAa> cos (<p - (l>)) 



GM+ 



\a(l -e 2 ) 



(Ae cos ((f> - to) + eAii) sin (<p - a>)) . 



In the limit e <sc 1 the relative velocity therefore becomes: 

Av 2 = v\ (Ae 2 + e 2 Aw 2 + S 2 AB ) . (26) 

The quantity Sab is defined within the coordinate system that 
is coplanar with orbit A. Since AO <K 1 and Aa « 1 we can 
apply the approximation introduced in Eq. (fTOt to express Sab 
as: 



S\ B = Ao- 2 + sin 2 (a)An 2 

and the relative velocity becomes: 

Av 2 = v\ [Ae 2 + e 2 Aco 2 + Aa 2 + sin 2 (a)AQ 2 ) . 



(27) 



(28) 



This result shows that relative velocities are generated for a 
variety of reasons. As in the coplanar case, relative velocities 
can be generated due to different eccentricities Ae, or phasing 
of pericentres Aa>. In three dimensions, relative velocities will 
also be caused by different inclinations Aa or phasing of nodal 
lines AQ. If Ae = Aa = 0, Au = AQ ~ 1 and sin (a) ~ a 



we recover the result from Lissauer & Stewart (1993) for ran- 
domised orbits: 



Av = vk Ve 2 + a 2 



(29) 



(25) 



4. Effect of disc gravity on planetesimal dynamics: 
single star case 

In this section we present the results of simulations of plan- 
etesimals that interact gravitationally with the disc and central 
star only. The planetesimals are on orbits that are inclined with 
respect to the disc midplane, and we switch off the drag force. 
The knowledge gained in this section will help to understand 
the results in later sections when the binary companion and the 
gas drag force are included in the model, and allow us to iso- 
late the effect that the disc gravity alone has on the dynamics 
of planetesimals on inclined orbits. 

Throughout this paper we use the figure convention that the 
top left panel is referred to as panel 1 , with the remaining panels 
being labelled as 2, 3, 4... when moving from left to right and 
from top to bottom. A simulation result for a reference particle 
is presented in FigQ](solid lines in panels 1-3 ). This planetes- 
imal has an initial semi-major axis of 10 AU, eccentricity of 
e -, =0.1 and relative inclination (inclination relative to the disc 
midplane) of S t =0.1. The quantities displayed are the nodal 
precession angle Q, (panel 1), inclination a, (panel 2) and apsi- 
dal precession angle (panel 3). The quantities are calculated 
with respect to a reference plane that is coplanar with the disc 
midplane (a, = Si). We observe that the disc gravity causes 
retrograde nodal precession (Q, < 0) about the disc angular 
momentum vector, and prograde apsidal precession (d>j > 0), 
while the inclination a, remains unaffected. 

We also studied planetesimals with different initial semi- 
major axes, eccentricities and relative inclinations. The results 
are summarised in panels 4-6 of Fig[TJ which show the apsidal 
(solid) and nodal (dashed) precession rates (a>j and Q,) as a 
function of: semi-major axis a,- (panel 4); eccentricity e,- (panel 
5); relative inclination S, (panel 6). The black lines show results 
for a disc model with the reference disc mass, and the red lines 
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are for a model with twice the reference mass. As expected, 
the precession frequencies (nodal and apsidal) scale roughly 
with the mass of the disc. Furthermore, the precession rates 
are larger in magnitude for particles that have smaller semi- 
major axes, as can be seen in panel 4. Panel 5 shows that the 
precession rates (nodal and apsidal) are weakly dependent on 
the eccentricity for the range in <?,■ shown. The dependence on 
relative inclination, however, is strong (see panel 6). 

So far we have discussed the evolution of orbital elements 
measured with respect to a reference plane that is coplanar with 
the disc midplane. Later on, however, we will consider plan- 
etesimal orbital elements measured with respect to the fixed 
binary plane, which will be inclined with respect to the disc 
midplane. If the inclination of the disc with respect to the bi- 
nary plane is smaller than the relative inclination of the parti- 
cles (ad < 6{), the nodal and apsidal precession rates caused by 
the disc remain unchanged, and panels 4-6 of Fig.Q~]apply. If the 
inclination of the disc with respect to the binary plane becomes 
larger than the relative inclination of the particles (ad > £,), 
however, the nodal and apsidal precession rates will be differ- 
ent. 

In order to understand the effect of the gravity of a disc 
that is highly inclined with respect to our reference plane, we 
transform the planetesimal position and velocity vectors via the 
transformation given by Eq. ©, using a representative inclina- 
tion of the disc of a d = y F = 0.78 (45°) > S t = 0.1 (5.7°) 
and assume that the disc is non-precessing with Q.d - Qft = 
(inclusion of a small precession frequency Qf + does not 
change the result.) The outcome of this transformation is shown 
by the dashed line in FigQ](panels 1-3) for the same reference 
particle as described earlier. We observe that the nodal preces- 
sion O, angle (panel 1) and inclination angle a, (panel 2) are 
oscillating around a fixed value. This can be understood as fol- 
lows. The angular momentum vector of the particle still pre- 
cesses around the angular momentum vector of the disc as be- 
fore. Unlike before, however, the disc has an inclination with 
respect to the reference plane of ad = 0.78 (45°) > 6j. Hence 
the precession of the particle angular momentum vector causes 
the measured inclination and nodal precession angles to oscil- 
late around those of the disc, with an amplitude given by the 
particle's relative inclination <5,-. We verify that for the refer- 
ence particle the nodal precession rate measured with respect 

to the disc midplane is Q 0.23 (-13.2°) per orbit, as seen 

from panel 4 (black dashed line). This corresponds to a preces- 
sion period of 27 orbits, coinciding with the observed period 
of oscillation of the inclination and nodal precession angles 
measured with respect to the reference plane (panels 1 and 2, 
dashed line). Additionally, the particle appears to precess apsi- 
dally in a retrograde sense (panel 3, dashed line). 

The precession rates (nodal and apsidal) are displayed as 
functions of semi-major axis, eccentricity and relative incli- 
nation in panels 7-9, respectively, for the highly inclined disc 
case, with the same line style and colour convention as used be- 
fore. As can be seen from these panels, the net nodal precession 
rate is zero, since the disc causes oscillation but no net change 
of the nodal precession angle. 

From panels 8-9 we observe that the apsidal precession rate 
is roughly independent of eccentricity e-, and relative inclina- 



tion Si. However it depends on the semi-major axis, a,, and also 
scales roughly with the disc mass, as shown in panel 7. For fu- 
ture purposes we fit this apsidal frequency by a second order 
polynomial: 



Md_ 
M 



C + C 



/ a 



\AU 
4.187- 10~ 2 



) + C2 (aTi 



i 

1Q- 3 and C 2 



(30) 



-4.242 • 



with Co = -4.187- 1Q- Z , Ci =5.118' 
10~ 4 . Pd is the orbital period at the disc outer edge (which is 
nominally located at r = 10 = 20 AU), Mj is the disc mass and 
M is the nominal disc mass for the reference model introduced 
in Sect. 3.3. The fit is displayed in panel 7 (dotted lines) and 
matches the numerical data (solid line) well. 

To summarise, the measured effect of the gravity of an in- 
clined disc on the planetesimal orbit depends on the ratio of 
the inclination of the disc with respect to the binary plane, aj, 
and the inclination of the particle with respect to the disc, 6,. If 
6i > ad the resulting nodal and apsidal precession rates caused 
by the disc are displayed in FigQ] (panels 4-6). If 5, < ad the 
precession around the disc angular momentum vector causes 
oscillations in the nodal precession and inclination angles with 
an amplitude that is given by The apsidal precession in this 
case is displayed in FigQ](panels 7-9) and can be approximated 
by Eq. 4H). 

5. Planetesimal dynamics in inclined binary 
systems 



Model 


7f 


M d 


D 


Si 


label 


[degrees] 


[M] 


[AU] 


[m] 


1 





1 


60 


10000 


2 


25 


1 


60 


10000 


3 


45 


1 


60 


10000 




45 


1 


60 


1000 




45 


1 


60 


100 


4 


45 


3 


60 


10000 


5 


45 


6 


60 


10000 


6 


45 


1 


90 


10000 


7 


45 


1 


120 


10000 


8 


25 


1 


60 


100 




25 


1 


60 


1000 




25 


1 


60 


10000 


9 


45 


1 


60 


100 




45 


1 


60 


1000 




45 


1 


60 


10000 



Table 2. Table of runs: the first column gives the run label; the 
second column gives the inclination of the binary companion 
to the gas-plus-planetesimal disc, jp ; the third column gives 
the disc mass in units of the reference mass; the fourth column 
gives the binary separation, D; the fifth column lists the plan- 
etesimal radii, Si, included in the models. 



We will now describe the dynamics of planetesimals that 
are orbiting in the full binary plus disc system. The model pa- 
rameters are summarised in table |2] The planetesimals are ini- 
tially set-up on circular, keplerian orbits that are coplanar with 
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Fig. 1. Upper panels: nodal (panel 1) and apsidal precession (panel 3) angles as well as inclination angles (panel 2) of a test 
particle experiencing the gravitational field of a disc and central star only. The solid lines show the results when viewed in a 
frame in which the disc is assumed to be coplanar with the reference plane. The dashed lines show the same quantities for the 
same model, but where the disc is now treated as if it was inclined by yp = 0.78 (45°) to the reference plane, as it will be 
when the binary companion is included. Middle panels: nodal (dashed lines) and apsidal (solid lines) precession rates when the 
disc is highly inclined as a function of semi-major axis a, (panel 4), eccentricity e,- (panel 5) and relative inclination 6/ (panel 
6) for a standard disc mass (black) and a disc with the double mass (red). Bottom panels: nodal (dashed) and apsidal (solid) 
precession rates caused by the highly inclined disc as a function of semi-major axis a, (panel 7), eccentricity e, (panel 8) and 
relative inclination <5, (panel 9) for a standard disc mass (black) and a disc with double that mass (red). In panel 7 a polynomial 
fit to the apsidal precession rates is shown (dotted lines). 



respect to the disc midplane (i.e. = 0). They are distributed 
radially in the interval a, e [6, 16] AU for models 1-7 (where 
the disc truncation radius is 20 AU). In an additional series 
of runs we confine the planetesimals to a narrow annulus with 
Aa = 10~ 3 AU centred around a = 10 AU (models 8 and 9) to 
maximize the number of collisions and study collisional veloc- 
ities in more detail. 

We treat model 3 as our reference case. The inclination of 
the binary companion to the disc is initially set to yp = 0.78 
(45°) and its separation from the central object is set to D — 
60 AU. The disc mass used in the reference model is = 
1M and corresponds to a scaled minimum mass solar nebula as 
explained in Sect. 3.3. In the other models we varied the initial 
inclination yp of the binary companion to the gas-planetesimal 
disc (models 1-3), the disc mass (models 4 and 5) and the 
binary separation D (models 6 and 7). 



5.1. Zero inclination case (model 1) 

In this model the binary orbital plane is coplanar with the 
planetesimal-plus-disc midplane. From time t — the plan- 
etesimals experience the gravity of the disc and gas drag forces 
(the planetesimals have size s, = 10 km), and the mass of the 
binary companion is increased to its final value over a time 
of four orbits. The evolution of the semi-major axes, a,, and 
eccentricities, <?,-, are shown in Fig|2] The colours correspond 
to different initial semi-major axes, with darker blue colours 
representing the inner planetesimals and the green-red colours 
representing planetesimals further out in the disc. We will also 
adopt this colour convention when discussing simulation re- 
sults later in the paper. As can be seen from Figj2] (panel 1), 
the binary companion does not change the semi-major axes of 
the bodies, as expected from secular theory. Eccentricities are 
generated, however, that are of the order of (e,- ~ 10~ 2 ), as can 
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Fig. 2. Orbital elements for the coplanar case. Panel 1 shows the 
semi-major axes of different planetesimals, where the colour 
convention introduced here is used in future plots. The eccen- 
tricity of the various planetesimals is shown in panel 2. The 
binary separaion D = 60 AU. 




20 40 60 
Time [orbits] 

inclinotion 




20 40 60 
Time [orbit] 

relotive inclinotion 




20 40 60 
Timeo [orbits] 



20 40 60 80 
Time [orbits] 



be seen from Figj2] (panel 2). Very modest eccentricities are 
generated because the planetesimals are orbiting in a slightly 
non-keplerian potential due to the disc gravity. But we also see 
that the eccentricities grow on a time scale of ~ 5 orbits due in- 
teraction with the binary whose mass grows over this time. The 
values of eccentricity ob t ained are consistent with those pre- 
sented bv lCiecieiag et al.1 d2007l) . who considered the evolution 
of planetesimals in a gas disc with a circular binary companion, 
but neglected the effects of disc gravity. Because of the closer 
proximity to the binary companion, the outer planetesimals are 
more strongly affected by these perturbations and their eccen- 
tricities are raised to larger values. Note that the data in this 
and the following figures has been smoothed over a temporal 
window of 1.8 orbits, but it is clear that after an initial rise, the 
eccentricities remain essentially constant for 80 orbits (7155 
years). 

Secular perturbations from an eccentric, coplanar binary 
companion lead to the generation of well-defined forced ec- 
centricities, and can also lead to strong alignment of perias- 
tra for same-sized plane te simals in the presence of gas drag 
dMarzari & Scholll l2000t iThebault. Marzari & Scholl \20o4 . 
This causes a dramatic reduction in the impact velocities for 
these planetesimals. But this effect is absent for a circular com- 
panion as considered here, and eccentricities are largely gen- 
erated by high frequency terms in the disturbing function. We 
find that the periastra are not well aligned near the beginning of 
the simulation, since a companion on a circular orbit cannot de- 
fine a direction of preferred alignment (this effect may be seen 
in Fig. 7 later in the paper where we plot results for a simula- 
tion with a low inclination (25°) binary compani on). This basic 



pint i s in agreement with results presented by ICiecielag et al 



poin 
(20C 



(2007), who also considered a circular binary. As such a binary 
companion on a circular orbit appears to be singular in its effect 
on the orbits of planetesimals embedded in a gas disc. This has 
potentially important consequences for the collisions between 
planetesimals reported later in this paper. 



Fig. 3. Orbital elements for the low inclination case jf — 25°. 
The colours represent planetesimals at different semi-major 
axes with the convention adopted from Figj2] The eccentricity 
is depicted in panel 1 . In panels 2 and 3 the nodal precession 
and inclination angles are shown, where the short dashed line 
represents the inner and outer edge of the disc. In panel 4 the 
relative inclination with respect to the disc is shown, where the 
short dashed line represents one pressure scale height and the 
dashed-dotted line indicates where 6, - aj. The binary separa- 
tion D = 60 AU. 



5.2. Low inclination case (model 2) 

In this section we describe the simulation results for model 2, 
for which only planetesimals of size 10 km were considered 
and the binary inclination yp = 25°. Figf3]shows the evolution 
of eccentricities <?,■ (panel 1), nodal precession angles O, (panel 
2), inclinations », (panel 3) and relative inclinations (5, (panel 
4) for the different planetesimals using the same convention 
of colours as described in Sect. 15.11 As in the coplanar case, 
eccentricities are raised by interactions with the binary com- 
panion, but because the binary orbit is circular the forced ec- 
centricity predicted by secular theory is negligible, and the ec- 
centricities are generated by high frequency perturbations. We 
can observe that these are somewhat lower than in the coplanar 
case due to a reduced coplanar component of the companion's 
gravity. 

In panel 2 of Figf3] we see that most of the planetesimals 
precess nodally at a joint rate with the disc, except for the out- 
ermost particle (red line), which will be discussed later in this 
section. This joint precession can be explained by the presence 
of the disc, with its gravitational field playing the major role. To 
illustrate this we also performed simulations of planetesimals 
that interact with the binary system only. The corresponding 
nodal precession and inclination angles of the planetesimals for 
such a simulation are shown in Figj4] Due to the secular per- 
turbation of the secondary star alone, the inclinations are ex- 
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pected to stay constant while the orbita l planes of the planetes- 
imals precess at their free particle rate (Pap aloizou & Terquem 
1995): 



free 



3 GMj 
80^ D 3 



(3 cos 2 a/ - 1). 



(31) 



which is a strong function of semi-major axis, as Qk ~ ^ - is 
the keplerian angular velocity. The left panel of Figj4] shows 
that the nodal lines of the planetesimals become progres- 
sively randomised, leading to the dispersion of the planetesimal 
disk into a cloud of bodies surro unding the star, as found by 
Marzari, Theb ault & ScholJ ( 2009 ). Comparing this behaviour 
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Fig. 4. Simulation results of planetesimals interacting with the 
binary system only (no gas disc). The left panel shows the nodal 
precession, and the right panel shows the inclination. 

with the results displayed in Fig|3](panels 2 and 3) we can see 
that the presence of the disc and its gravitational field causes 
fundamentally different behaviour of the planetesimals. The 
planetesimals remain coupled to the disc, as their nodal pre- 
cession angles stay close to the disc precession angle (dashed 
line in panel 2 of Fig. [3}. The role of gas drag in determining 
the relative inclination between the disc midplane and the plan- 
etesimals is negligible for this run for which the planetesimal 
size is 10 km. But for smaller planetesimals gas drag does play 
a role, as we will describe later in the paper. 

In FigOwe show the time evolution of the nodal preces- 
sion rates for the particle at r — 13.3 AU (panel 1) and r = 6.2 
AU (panel 2). The various curves represent contributions due 
to gas drag (blue line), disc gravity (red line), binary compan- 
ion (green line) and the total rate (black line) as calculated in 
Sect. 3.4. Because of its different radial location, the outer par- 
ticle orbit precesses faster in the retrograde sense than the inner 
particle if the disc is absent. Confirming this we can observe in 
Fig[3(panels 1 and 2) that Q, due to the companion (green line) 
is about O, = -0.09 (-5.15°) per orbit for the outer body, and 
CI, - -0.03 (-1.7°) per orbit for the inner body. It is the disc 
gravity, however, that compensates for this effect and causes 
the planetesimals to precess at a common net rate. Since the 
disc precesses rigidly at a rate &d Qf = -0.06 (-3.4°) 
per orbit, which corresponds to the free particle rate at about 
cii — 11 .5 AU, disc gravity (red line) gives a positive contribu- 
tion to the precession rate for the outer particle (Q,- > 0) and a 
net negative contribution for the inner particle (£% < 0). Since 
the particle-plus-disc system precesses in the retrograde sense 



the outer particle precession is slowed down, while the inner 
particle precession is speeded up by the disc, such that both 
bodies precess together with the disc on average at the rate Clf. 
This result illustrates the fundamental importance of including 
the effects of disc gravity when calculating the dynamics of 
planetsimals in misaligned binary systems. 

The nodal precession and inclination angles of the planetes- 
imals show oscillations, as seen in Figj3](panels 2 and 3). These 
are caused by the disc gravity, as explained in Sect. 4 and illus- 
trated by the red line in panel 5 of Fig. [5] The planetesimals 
effectively precess around the disc angular momentum vector 
at the same time as the disc-plus-planetesimal system precesses 
around the binary angular momentum vector. Initially the plan- 
etesimals are coplanar with the disc midplane ((5,- = 0), but as 
the system evolves during early times, differential nodal preces- 
sion between particle i and the disc occurs, leading to a build 
up of relative inclination (5,, since Sj oc (Q,- - Q d ) 2 as shown 
by Eq. (TTOb . The differential precession is greatest for planetes- 
imals whose semi-major axes deviate most from a = 1 1 .5 AU 
(the radius at which planetesimals naturally want to precess at 
the same rate as the disc), so the relative inclination is largest 
for particles that are furthest inside (black line) or outside (yel- 
low line) this value, as demonstrated in panel 4 of Fig. [3] 

The disc model has an aspect rato H/R = 0.05, so that 
planetsimals which develop relative inclinations > 0.05 will 
spend large fractions of their orbits away from the disc mid- 
plane. Planetesimals interior to 8 AU and exterior to 13 AU 
have relative inclinations <5; 0.1, and so spend large fractions 
of their orbits essentially outside of the disc. We also observe 
that the oscillation periods of the inclination, c,-, and nodal pre- 
cession angles, Q,-, varies among the planetesimals. This is ex- 
pected due to their different semi-major axes and relative incli- 
nations, and the observed periods are all consistent with FigQ] 
(dashed lines in panels 4-6), which we recall shows how the 
disc alone acts on the planetesimals. 

We note that the inclination oscillations, seen in Fig|3] 
(panel 3), begin in opposite senses for planetesimals whose 
semi-major axes are above or below a = 1 1 .5 AU (i.e. plan- 
etesimals below a - 1 1 .5 AU are perturbed onto higher incli- 
nation orbits, while planetesimals exterior to a = 1 1 .5 AU are 
perturbed onto lower inclination orbits). Panel 3 of FigJ5]plots 
the inclination angle difference a, - ad versus the precession 
angle difference Q, - for a planetesimal located at 13.3 AU. 
A similar plot for a planetesimal at 6.2 AU is shown in panel 4. 
Each plot shows the trajectory of the tip of the planetesimal an- 
gular momentum vector relative to the disc angular momentum 
vector (which is located at the origin). For the inner particle 
the companion induces a differential precession Q,- - > 
and the anti-clockwise precession around the disc angular mo- 
mentum vector causes the planetesimal to approach a higher 
inclination orbit a, - > during the first half of this preces- 
sion cycle. For the outer planetesimal the induced differential 
precession is Q,- - Q,/ < 0, and the anti-clockwise precession 
causes perturbation onto a lower inclination orbit a, - < 0. 

We now return to the planetesimal that is orbiting at a — 15 
AU. This body shows quite distinct behaviour from all the other 
bodies, as shown in Figj3](red lines). This particle is dominated 
by the gravity of the companion, and it nodally precesses fast 
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Fig. 5. Panels 1 and 2 show the nodal precession rates of a plan- 
etesimal with semi-major axis at a, = 13.3 AU (panel 1) and 
a, = 6.2 AU (panel 2). Shown are the contributions of the disc 
gravity (red line), binary companion (green line), gas drag (blue 
line) and total rate (black line). In panels 5 and 6 we also show 
the rates of inclination change with respect to the binary and 
disc midplane respectively for the inner particle. In panels 3 
and 4 the trajectory of the tip of the particle's angular momen- 
tum vector is shown projected onto the (£2, - D.^, a-, - a,j) plane 
for the outer (panel 3) and inner particle (panel 4). 



enough to decouple from the disc. As it precesses away from 
the disc, the relative inclination Si grows (panel 4 of Fig|3]l. The 
precession rate around the disc angular momentum vector is a 
decreasing function of relative inclination, so the precession 
period becomes long for this particle. Because Q, - £2</ < 0, we 
observe the quasi-monotonic decrease of inclination in panel 
3 of Figj3] until the reversal point at t = 47 orbits, when the 
differential nodal precession is £2, - = n. The planetesi- 
mal orbital plane then approaches the disc midplane from the 
other side (Q, - > 0), and we observe the quasi-monotonic 
increase of inclination after t = 47 orbits. 

Over longer time scales than we have been able to consider 
because of the computational expense of running the simula- 
tions, we expect the outermost planetesimals which decouple 
from the disc to show continued oscillations in their inclina- 
tions, with an oscillation period equal to the synodic preces- 



sion period. Ignoring the possible growth of planetesimals into 
planets, or their collisional destruction, those bodies which re- 
main almost coplanar with the disc should continue to do so 
over long time scales until the disc mass decreases due to vis- 
cous evolution and/or photoevaporation. Significant reduction 
of the disc mass would eventually allow the planetesimal dy- 
namics to become dominated by the companion star, and they 
would decouple from the disc. But, we also note that that the 
disc and binary orbit plane will evolve toward alignment on 
the viscous evolution time at the outer edge of the of the disc 
dTerquem et al.lll999t ILarwood et al.|[l99oi iFragner & Nelson 



2010). So it is possible that an initially misaligned protoplane- 
tary disc may align significantly over its lifetime, bringing with 
it any planetary system which has formed within it. The final 
degree of misalignment for a planetary system formed in an in- 
clined protoplanetary disc may therefore be substantially less 
than the initial misalignment of the protostellar disc. 

5.2.1. Collisional velocities in low inclination case 
(model 7) 

In this section we present the results of a simulation which ex- 
amines collisions between planetesimals in a system where the 
binary companion has an inclination of yp = 25° relative to 
the disc midplane. To increase the collision rate to statistically 
meaningful values which can be measured in a simulation, we 
initialise the planetesimals with a narrow range of semi-major 
axes (Aa = 10 ^ AU) centred around a = 10 AU. We consider 
three different planetesimal sizes (100m, 1km and 10km), and 
for each size we include 50 particles. We check the condition 
for orbital crossing given by Eq. ([18} every 100 time steps for 
each of the 11175 particles pairs (100 time steps corresponds 
to 0.01 1 orbital periods at 10 AU, the orbital radius of the plan- 
etesimals). If the orbit crossing condition is satisfied, then we 
calculate the velocity at the crossing point for both planetesi- 
mal orbits according to Eq. ([T9V 

Before discussing the results of model 7, it is worth recap- 
ping what we might expect based on previous work in which 
the gravity of the disc was neglected. Same-sized planetesi- 
mals being perturbed by an eccentric, coplanar binary com- 
panion will experience a growth in their forced eccentricity, 
but gas drag damping will cause stron g orbital phasing dramat- 
ically reducing collisional velocit ies (IMarzari & Scholl 2000; 
iThebault. Marzari & Scholll 120061) . The different phasing of 
pericentres for planetesimals of different size, however, leads to 
large collision velocities which are likely to be disruptive. The 
inclusion of a small inclination (a, < 5°) causes different sized 
planetesimals to orbit in different planes, such that collisions 
between similar sized bodies are more frequent than between 
different sizes. The fact that the pericentre phasing is main- 
tained in this scenario means that planet esimal growth ma y be 
more likely to occur in inclined systems (Xie & Zhoull2009l) . 

The collision velocities we obtain are shown in Fig|6](solid 
lines) for collisions between equally sized (panel 1) and differ- 
ently sized planetesimals (panel 2). As mentioned in Sect. 15. II 
an important difference between our set-up and previous work 
is that we consider a circular, inclined binary companion, 
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Fig. 6. Average collisional velocities in the low inclination yp = 25° case between equally sized (left panel) and differently sized 
planetesimals (right panel) for planetesimals centred around 10 AU from the central star. Left panel: collisions between 10km- 
10km bodies (green-solid); lkm-lkm (blue-solid); lOOm-lOOm (red-solid). Right panel: collisions between lOkm-lkm bodies 
(green-solid); lOkm-lOOm (blue-solid); lkm-lOOm (red-solid). Threshold velocities for catastrophic disruption corresponding to 
the different size-combinations for the weak (dash-dotted line) and strong aggregates (dashed line) are also shown. 



resulting in a broad distribution of planetesimal longitudes 
of pericentre, w,-, even for planetesimals of the same size. 
Consequently, we see that the collision velocities for equal 
sized bodies in panel 1 are quite large, being between 50 - 
70 ms"' at the end of the simulation. Although it appears that 
the circular binary is largely responsible for the growth of ec- 
centricity and the misaligned periastra of orbits for same sized 
planetesimals, it is possible that gravitational perturbations as- 
sociated with spiral waves in the disc also contribute. The col- 
lision velocities for differently sized bodies, however, are even 
larger than for same-sized bodies, and exceed 400 ms _1 after 
80 orbits (see panel 2). 

Fig(7] shows averages of all the quantities that determine 
the relative velocities according to the analytical estimate given 
by Eq. d28l i. and we can use this data to determine the main 
contributors to the collision velocities shown in Figj6] Panel 1 
shows the difference of longitude of pericentre Acjj, panel 2 
shows the average eccentricity <?, panel 3 shows the difference 
in eccentricity Ae,, panel 4 shows the difference in inclination 
Att;, panel 5 shows the angle between the nodal lines A£2, and 
panel 6 shows the average inclination <?,. Note that these aver- 
ages were obtained only for orbits which mutually cross, and 
do not represent the distributions for the whole ensemble of 
particles. Using the numbers extracted from these figures in 
Eq. ( l28b we can reproduce the relative velocities shown Figj6] 
with good accuracy, indicating that Eq. (|28| | is a valid approxi- 
mation and that our collision detection technique is generating 
collision velocities which are consistent with expectations. 

For collisions between planetesimals of the same size the 
dominating contribution to the square of the relative velocity 
comes from the term which is proportional to e^Au^, while all 



other terms are at least an order of magnitude smaller. This im- 
plies that planetesimals with the same size all orbit more or 
less in the same plane (as shown by panels 4 and 6), and rela- 
tive velocities are generated due to misalignment of their orbit 
pericentres, Aw,, a result which is broadly consistent with those 
obtained bv lXie & Zhoul (120091) . We note that our collision de- 
tection method is not strictly valid in this regime, because it 
predicts the locations of the orbit crossing points incorrectly 
for coplanar orbits. But, because the average pericentre mis- 
alignment is on the order of unity (measured in radians), the 
average collision velocities reported by the algorithm are con- 
sistent with the analytical prediction v co u - e Aw to within a 
factor of order unity. Under these conditions, incorrect detec- 
tion of the orbit crossing locations still leads to estimates of the 
average collision velocities being approximately correct. 

Although panel 1 of Fig|7] shows that the misalignment of 
pericentres increases slightly over time when considering only 
those particles whose orbits intersect, we find that the average 
pericentre misalignment calculated over all same-sized particle 
pairs actually decreases during the simulation. The reason for 
this is not clear, but it is interesting to note that even when 
the binary orbit is circular, the long term evolution appears to 
be toward one where the pericentres tend toward alignment, 
and this effect is most marked for the smallest planetesimals 
suggesting that it is an effect due to gas drag. 

As seen in panel 2 of Fig[6j planetesimals of different sizes 
tend to have much larger collisional velocities. One reason 
is the increased misalignment of pericentres, as can be ob- 
served in Fig[7] (panel 1, dashed lines). This may be a di- 
rect consequence of gas drag, as discussed in previous work 
dMarzari & Scholll l2000l iThebault. Marzari & ScholU |2006|) . 
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Fig. 7. Average orbital parameters of the colliding planetes- 
imal pairs in the low inclination yp = 25° case. The line 
colours and styles are as follows: green-solid (lOkm-lOkm); 
blue-solid (lkm-lkm); red-solid (lOOm-lOOm); green-dashed 
(lOkm-lkm); blue-dashed (10km- 100m); red-dashed (1km- 
100m). 



but it may also be due to the fact that particles which expe- 
rience stronger gas drag orbit closer to the disc midplane than 
larger planetesimals do. This then affects the gravitational per- 
turbations experienced by the particles as they travel through 
the disc, which may increase the pericentre misalignment. But 
a more important contribution to the square of the relative ve- 
locities comes from the term oc sin 2 (a,) AO 2 . We find that due 
to the different gas drag strengths, planetesimals of different 
sizes tend to precess nodally at different rates for some initial 
time until disc gravity causes them to precess together with the 
disc. While smaller planetesimals tend to precess together with 
the disc immediately and consequently their relative inclina- 
tions stay low, larger planetesimals tend to precess at their own 
free particle rates for a longer time, causing a larger build up 
of relative inclination. In other words, although the planetesi- 
mal swarm as a whole is forced on average to precess with the 
disc by the disc gravity, smaller planetesimals oscillate about 
the midplane of the disc with a smaller amplitude than larger 
planetesimals. This induces a large misalignment of their or- 



bital planes A£2, as seen in panel 5 of Figj7] (dashed lines) in 
addition to the misalignment of their pericentres. 

The frequency of collision between bodies is an impor- 
tant fa ctor during the growth of planetesimals, and lXie & Zhoul 
(2009) lave suggested that the differential phasing of nodal 
lines for planetesimals of different sizes may increase the rel- 
ative importance of collisions between similar sized objects 
rather than different sized bodies. We have examined the colli- 
sion frequency reported by our collision detection technique 
by simply counting the number of orbit crossing events re- 
port ed during the simul ation. In basic qualitative agreement 
with IXie & Zhoul ([2009), we find that the collision frequency 
during the simulation between same-sized objects is a factor of 
3-6 times larger than between different sized objects. Our re- 
sults thus support the idea that planetesimal growth in inclined 
binaries is likely to proceed via accretion of similar sized bod- 
ies. 

We conclude that relative velocities between differently 
sized planetesimals tends to be large in binary systems whose 
orbital plane is misaligned with the planetesimal-plus-disc 
plane, while relative velocities between same sized planetesi- 
mals are largely unaffected by this non-coplanarity. But the cir- 
cular orbit of the binary companion prevents strong pericentre 
alignment for similar sized bodies, leading to significant col- 
lisional velocities in this case too, in contrast to the situation 
observed when the companion is on a coplanar eccentric orbit. 
The question of what happens in the case of an eccentric, in- 
clined companion unfortunately goes beyond the scope of this 
paper but should obviously be the focus of future work. 

We now want to determine whether the collisional veloc- 
ities seen in Fig|6]will lead to growth or erosion of the plan- 
etesimals. In principle there are three possible outcomes: accre- 
tion, in which the largest body involved in the collision contains 
more mass after the collision; catastrophic disruption, in which 
more than half of the total mass of the system is dispersed af- 
ter the collision; erosion (or cratering), where the largest body 
loses a small amount of mass during the collision. Only accre- 
tion leads to the growth of a planetesimal. 

For simplicity in interpreting our results, we adopt 
the universal law for the largest remnant mass from 



Stewart & Leinhardt (2009): 



Mir 

M, nt 



2Q D ' 



(32) 



where M\ r is the mass of the largest post-collision remnant, and 
M, , - M\ +M2 is the total mass of the two colliding bodies Mi 
and M2. The quantity Qr = 0.5 ^^- Av 2 is the reduced mass 
kinetic energy scaled by the total mass of the colliding system. 
For equally sized bodies accretional collisions require > 
5, from which it follows that Qr < Qo- Hence Qd is called 
the catastrophic disruption limit of the collisions (the energy 
required to disperse half the total mass). If Qr > Qo collisions 
between equally sized bodies lead to catastrophic disruption. 
When considering collisions between differently sized bodies 
the condition has to be modified. Let Mi » M2 such that M2 = 
//Mi with n <K 1 . The condition for accretion is now M/ r > Mi , 
from which it follows that Qr < Qd, which implies for 
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Fig. 8. Threshold velocities for collisions between equally 
sized planetesimals as a function of their spherical radius R\, 
The disruption curves for the weak (dashed-dotted line) and 
strong rock (dashed line) are shown for two different impact 
velocities Vi = 50m/s (green) and Vj = 2km/s (blue) along 
with the solution for rubble piles (black-solid line). The bodies 
escape velocity (black-dotted line) is also shown. 



the relative velocity threshold: 



6v D = 2 vrr^ Ve^ 



(33) 



The disruption limit Qd is sensitive to factors that influence the 
energy and momentum coupling between the colliding bodies 
(i.e. impact velocity, m ass ratio and material proper ties such as 
strength and porosity'). IStewart & Lein hardt (2009) show that 
the catas trophic disruption curve Qp can be fit to an analytical 
formula dHoussen & Holsappld[l990. 19991) : 



Qd - ylsK-12 + ?C K i2 



<]vf 



(34) 



where fiM and 4>m are material properties. In this expression R\2 
is the spherical radius of the combined mass M m assuming a 
density of p, = 1 g cm 4 . Since we use a density of p s = 2 
g cirT 3 this radius is given by R12 — 2^(1 + fi)*Ri, where Ri 
is the spherical radius of the larger mass Mi involved in the 
collision. The first term on the right hand side of Eq. (l34l rep- 
resents the strength regime, while the second term represents 
the gravity regime. Very small bodies are supported by mate- 
rial strength (first term), which decreases as the planetesimal 
size increases. As gravity becomes more important for larger 
bodies the second term becomes dominant and increases the 



disruption limit again. IStewart & Leinhardtl (120091) derive the 
constants qs , qc,Hiw and 4>m for weak aggregates, such as weak 
rock and porous glass, by fitting the disruption curve to their 
numerical results. For weak aggregates they find (in cgs units) 
qs = 500, qc = 10~ 4 , hm - 0.4 and (p M = 7. For strong rocks 
they use the basalt laboratory data and modelling results from 



iBenz & As nhaug 
jUm = 0.5 and <pM = 8. In the limit of very large bodies the 
disruption curve can be best fitted by their resul t s of colliding 
rubble piles. In this regime IStewart & Leinha rdt (2009) find 



Qrp = 1.7-10 



6r>2 

K 12 



(35) 



for equal mass bodies and a factor of about three larger than this 
if the impactor size is much smaller than the target size. For il- 
lustrative purposes we present the relative velocity threshold 
6vp obtained by combining Eq. d33l and Eq. (t34b in the case 
of identical colliding bodies (fi = 1) in Fig|8] Note that in the 
limit of large bodies at low impact velocities, we set Qp = Qrp, 
because the rubble-pile solution for equal mass bodies defines 
a lower limit for disruption in this regime. The figure shows 
the limiting velocity below which collisions will lead to net 
accretion as a function of the planetesimal radius Ri for the 
weak (dashed-dotted line) and strong aggregates (dashed line) 
at two different impact velocities Vj = 2km/s (blue line) and 
V/ = 50m/s (green line). It is apparent that the threshold ve- 
locity is increased for larger impact velocities as more energy 
is partitioned into shock deformation. In contrast, at low im- 
pact velocities the momentum coupling between the colliding 
bodies is more efficient and less energy is needed to cause ero- 
sion of planetesimals. As gravity becomes more important than 
strength for very large bodies the curves join onto the disrup- 
tion curve for colliding rubble piles (black solid line). In this 
regime the threshold velocity is about a factor of 3-5 larger 
than the escape velocity from the bodies surface (dotted line). 
To compare the relative velocities obtained from our numerical 
simulations with the disruption threshold velocity, we evaluate 
Eq. (f33j) for equally (ji — 1) and differently sized fj. = 10~ 3 
planetesimals using V/ = Av in Eq. (l34l and choosing Ri as the 
bigger of the two planetesimals involved in the collision. If the 
rubble-pile disruption limit gives a higher estimate for a given 
size combination we use this limit instead. 



We note that lStewart & Leinh ardt (2009) considered target- 
projectile mass ratios only down to 0.03, instead of the val- 
ues 10 3 and 10 6 that apply to collisions of 1km or 100m 
sized b odies with a 10km sized tar get. As such we are apply- 
ing the IStewart & Leinhardtl (120091) results in a regime where 
catastrophic disruption is unlikely to occur, and therefore out- 
side of the regime of validity of their study, strictly speaking. 
Nonetheless, it reasonable to assume that even if catastrophic 
disruption does not occur during high velocity impacts, accre- 
tional growth is also unlikely due to erosion. 

The results are plotted in Figj6]for the strong (dashed line) 
and weak aggregates (dashed-dotted line). Except for colli- 
sions between equally sized 10km bodies (green), the colli- 
sional velocities are always substantially larger than the ero- 
sion/disruption threshold. Hence we conclude that growth of 
planetesimals probably can not occur for planetesimal sizes 
< 10km. For the equally (10km) sized bodies on the other 
hand the situation is not as clear. Although the collisional ve- 
locities lie below the strong aggregate threshold (green dashed 
line) they are also slightly above the weak aggregate thresh- 
old (dashed-dotted line). This suggests that collisions might 
also lead to erosion in this size regime. But it should be noted 
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that we have not included a realistic planetesimal size distri- 
bution in these simulations, and it is possible that collisions 
between planetesimals clustered in size around 10km, whose 
orbit planes are very modestly inclined with respect to one 
another, may allow accretion to occur for strong aggregates. 
Examination of this possibility will require future simulations 
that adopt a more realistic and continuous size distribution. 

If collisions between 10km sized planetesimals do lead to 
growth then this will most definitely not occur in the standard 
runaway regime, because the escape velocity is ~ lOm/s. It 
is possible, however, tha t type II runaway growth may occur 
(Kortenkamp et al. 2001) in which the large velocity dispersion 
causes growth to be orderly while the bodies remain small, but 
enters a runaway phase when large bodies form whose escape 
velocities are larger that the velocity dispersion. We conclude 
that in an inclined binary system on a circular orbit, relative 
velocities are excited that will lead to erosion of planetesimals 
for sizes < 10km which we have considered here, and growth 
is uncertain for 10km sized bodies. 

5.3. High inclination cases and the Kozai effect 

For large inclinations between the binary and planetesimal or- 
bital planes, it is possible that the Kozai effect will switch on, 
causing cy clic variation s of planetesimal eccentricities and in- 
clinations ( Kozailll962h . In this section we describe this effect 
in more detail in order to simplify the understanding of results 
which are described later in this paper. The equations describ- 
ing the secular evolution of the inclination a,-, e ccentricity, e,, 
and lo ngitude of pericentre, w,, can be written as (Inn anen et al 
1997b : 

15 sin (2<y,) sin (a,) cos (ar,) 



dost 

7fr 



(36) 



1 - e? sin (2o>i) sin (a,-) 



(37) 



de, 15 f 

^7 = ^VW{2(l-^) + 5sin 2 ( W ,)[^-: 

where for simplicity of notation we have introduced the time 
variable: 



1 («;)]} (38) 



t = [D^KilGMeYh. 



(39) 



Hence the time scale on which any of the given quantities 
change is oc D 3 . It can be seen from Eq. (1381 that there is a 
value of ojj for which <i», = 0. To lowest order in e,- this occurs 
when 

2 



sin (u>i) 



(40) 



5 sin (a,o) 

where a,o is the initial inclination of planetesimal i. Having ob- 
tained this value, the apsidal precession is halted, as w, = 0. 
This leads to the exponential growth of eccentricity. To illus- 
trate this we can substitute Eq. (1401 into Eq. d37l > to find: 



dei 
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— = —sin (a i0 ) sin (w ; ) cos (w,)e,- 
or 4 




sin (ff, ) - - 



(41) 



Hence the Kozai effect is only active for initial inclinations 

sin (atjo) ^ -^f or written differently a,o ^ 0.68 radians. The 
critical angle for which the Kozai effect switches on is thus 
ok = 39.2°. Additionally it can be shown that the Delaunay 
variable, D,, (which is equivalent to the component of the plan- 
etesimal angular momentum which is parallel to the binary an- 
gular momentum vector) is a constant of the motion, where D, 
is defined by 



(1 - ef) cos (or,-)- 



(42) 



As secular variations do not change the semi-major axes, this 
implies that an increase in eccentricity is coupled to a decrease 
in inclination, and vice versa. 

To study the Kozai effect in the absence of the disc, we 
have performed N-body simulations of planetesimals inter- 
acting with the binary system only. The results are depicted 
in Fig|9]for planetesimals with various initial inclinations. It 
can be seen that the eccentricities and inclinations undergo 
cyclic variations for planetesimals whose initial inclination are 
a-, > ax (yellow, green and light blue lines in Figj9j. This can 
be explained as follows. As the eccentricity grows and the in- 
clination decreases (panels 1 and 3), the latter eventually hits 
the Kozai threshold, ckj^, at which point there is no value for 
for which expression (1401 holds. Consequently w,- + and co, 
evolves to values for which <?,■ < and a-, > 0. Hence, even- 
tually a>i = can be obtained again, but now for a value of 
which is causing the opposite evolution of eccentricities and in- 
clinations until the original configuration is recovered and the 
Kozai-cycle is completed. We can clearly see in Fig(9]how the 
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Fig. 9. Orbital elements of planetesimals with different initial 
inclinations interacting with the binary system only. Panel 1 
shows the eccentricity and panel 2 shows the nodal precession 
angle. The inclination with respect to the binary plane and the 
longitude of pericentre (apsidal precession) are depicted in pan- 
els 3 and 4, respectively. 
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increase or decrease of eccentricities and inclinations is con- 
nected to the halting of apsidal precession during each half- 
cycle of the Kozai mechanism (panel 4). Additionally, we ob- 
serve that the time scale on which the eccentricities and incli- 
nations vary are increased as the initial inclination a® is de- 
creased. Once the initial inclination is a® < ax (Fig|9j blue 
and black lines) the condition <w, = can never be obtained. 
Consequently there are no net changes in eccentricities and in- 
clinations and the Kozai effect is switched off. 



5.4. High inclination case (model 3) 

In model 3 we increased the inclination between the disc- 
planetesimal system and the binary plane to yp = 0.78 (45°). 
As described in the previous section, once the inclination ex- 
ceeds a critical value ok = 0.68 (39.2°), the Kozai effect may 
lead to large eccentricity and inclination variations of the plan- 
etesimals, provided the disc mass is not too large (see later for 
a discussion about the conditions under which the Kozai ef- 
fect operates). The presence of our scaled minimum mass solar 
nebula disc does not prevent the onset of the Kozai effect for 
planetesimals with semi-major axes a < 12 AU, as seen from 
FigJTU] (panels 1 and 3), where we use the same colour con- 
vention for depicting the planetesimals as a function of initial 
semi-major axis as in models 1 and 2. 
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Fig. 10. Orbital elements for the highly inclined case jp = 
0.78 (45°) (model 3). Panel 1 shows the eccentricity grows due 
to the Kozai effect for all but the outermost planetesimals . The 
nodal precession and inclination angles are depicted in panel 2 
and 3, respectively, where the short dashed lines represent the 
inner and outer disc edge. The long dashed line in panel 3 rep- 
resents the threshold inclination of - 39.2° above which 
the Kozai mechanism is expected to operate. The apsidal pre- 
cession angles are depicted in panel 4. 



The eccentricity grows earliest for planetesimals with 
larger semi-major axes (panel 1, green lines). When the Kozai 
mechanism starts operating, the inclination starts to decrease. 
In the absence of the disc, this decrease would continue until 
the threshold of ax = 0.68 (39.2°) is approached, after which 
the inclination would increase and the eccentricity would de- 
crease again. In the presence of the disc, however, the situation 
is different. As the Kozai effect reduces the inclination, a, the 
nodal precession induced by the binary accelerates according 
to Eq. (t3TT >. This causes the planetesimal orbits to become sig- 
nificantly inclined relative to the disc midplane (see panel 3 of 
FigflOb. Concurrent precession about the disc angular momen- 
tum vector due to the disc gravity causes a quasi-monotonic 
decrease of inclination relative to the binary. The planetesi- 
mals are thus perturbed by the disc onto orbits with inclination 
a, < ax, causing the Kozai effect to switch off. At this stage 
the planetesimal has only completed a portion of its Kozai cy- 
cle and is left with high eccentricity (panel 1) and a reduced 
inclination (panel 3), at least for the duration of the simulation. 
The apsidal precession angles are depicted in panel 4 of FigJTOl 
During the time when the Kozai effect is operating (a, > ax), 
the apsidal precession is halted such that d>, = 0. After the 
Kozai effect has switched off (<?, < a^) we observe prograde 
apsidal precession (w,- > 0). Although this precession rate is 
dominated by the perturber the disc enhances the prograde pre- 
cession rate, since 5, > at this stage. The total rate is thus 
given by the sum of contributions by the binary companion and 
disc. 

Planetesimals with semi-major axes a, > 12 AU decou- 
ple from the disc due to the differential precession induced by 
the binary companion. Consequently their inclinations get per- 
turbed below (*k before the Kozai effect can start to operate, 
and their eccentricities stay low, at least for the duration of the 
simulations that we present here. 

The simulation presented in Fig[l0]was only run for 110 
orbits measured at the disc outer edge (~ 10 4 years). Those 
planetesimals which have experienced the Kozai effect have 
only undergone half a Kozai cycle, which stalls because the 
inclination falls below a^. But we see in panel 3 that the in- 
clination relative to the binary is increasing toward at the 
point when the simulation ends, because the planetesimals have 
precessed nodally by more than 180° with respect to the disc. 
We would therefore expect that the Kozai effect will switch 
on again when c,- > a^, allowing the previous Kozai cycle to 
complete. Similarly, we see that the inclinations of the plan- 
etesimals located beyond 12 AU are also increasing toward a^, 
such that the Kozai effect will switch on for them eventually. 
A long term effect of the disc thus appears to be to prolong the 
period associated with the Kozai cycles because of its effect 
on the planetesimal inclinations. An important consequence of 
this is that it increases the dephasing of the Kozai cycles at dif- 
ferent locations in the disc, ensuring that over long times plan- 
etesimals with different semi-major axes are at very different 
phases of their Kozai cycles. There will therefore be very large 
variations in both eccentricity and inclination within the plan- 
etesimal swarm, leading to very large collisional velocities be- 
tween planetesimals that are well separated in their semi-major 
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5.4.1 . Varying planetesimal sizes (model 3) 

In the previous discussion we focused on planetesimals with 
size si = 10km, for which the gas drag forces are very weak. 
As we decrease the size of the bodies in the system, the gas 
drag becomes more important since the stopping time is oc s,. 
Results for the 1km and 100m sized planetesimals are shown in 
the left and right panels of Fig[TT]respectively. The Kozai effect 
leads to the growth of eccentricities (panels 3 and 4) and rela- 
tive inclinations (panels 7 and 8) for both particle sizes shown. 
This is not surprising, since the time scale on which we expect 
substantial damping of eccentricity and inclination to occur is 
of the order of the gas drag stopping time. For our disc model 
this is given by: 



t s =0.81 



a Y 
AU/ 



Pd, 



(43) 



where s, is expressed in metres. For the 100m sized planetes- 
imals at a = 6-15 AU the stopping time is of the order of 
10 3 - 10 4 orbits, respectively, and is a factor of 10 longer for 
the 1km sized planetesimals. This is much longer than the time 
scale on which the Kozai effect operates (~ 10 2 orbits), and 
hence we would not expect the gas drag to prevent growth of 
eccentricity and relative inclination in this size range. We com- 
ment here (without showing results), that we also found the 
Kozai effect to operate for 10m sized bodies, for which the drag 
and Kozai time scales are similar. 

The increased effect of gas drag causes the semi-major axes 
to decay (see panels 1 and 2 in FigfTTb. As planetesimals ex- 
perience the Kozai effect, and develop large eccentricities and 
relative inclinations, their relative velocities with respect to the 
gas disc becomes very large (since |v; - v| 2 - v 2 K (e? + <5 2 ) from 



expression H29I ). As they travel through the disc they lose ki- 
netic energy rapidly to the gas, and their semi-major axes de- 
cay. We note that those 100m sized bodies lying beyond 12 AU, 
which do not experience the Kozai effect during the simulation, 
also undergo rapid decay in their semi-major axes. These par- 
ticles develop large inclinations relative to the disc due to the 
rapid nodal precession induced by the companion, and the re- 
sulting gas drag as they pass through the disc leads to their 
orbital decay. 

Before the Kozai effect starts to operate we can see that the 
increased effect of gas drag for the 100m sized bodies (right 
panels of Fig. ITTb causes damping of the relative inclination 
(panel 8) and therefore a decrease of the amplitude of the oscil- 
lation about the disc midplane caused by disc gravity (panel 6). 
In other words, the increased gas drag causes the 100m sized 
bodies to remain closer to the disc midplane, with the oscilla- 
tions in the amplitude of relative inclination being damped to 
a larger degree than occurs for larger bodies. This effect will 
become important when discussing relative velocities between 
the differently sized planetesimals, as the size dependence of 
the inclination perturbations will have an effect on the Kozai 
effect operating on the planetesimals, causing planetesimals of 
different sizes to orbit in planes which are significantly inclined 
with respect to one another. We note that the relative inclina- 
tions of both the 1km and 10km sized bodies are found to be 
similar, since the gas drag is weak in both of these cases. 



5.4.2. Collisional velocities for the highly inclined case 
(model 9) 

In this section we investigate the relative velocities for the high 
inclination (jf = 0.78) (45°) case. The procedure to detect 
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Fig. 11. Orbital elements for the 1km (left panels) and 100m 
sized planetesimals (right panels). The semi-major axes are 
shown in panels 1 and 2. In panels 3 and 4 eccentricity growth 
due to the Kozai effect is still apparent for both planetesi- 
mal sizes. The inclination is shown in panels 5 and 6, where 
the short dashed lines represent the disc inner and outer edge, 
and the long dashed line marks the threshold inclination above 
which the Kozai effect operates. In panels 7 and 8 the relative 
inclination with respect to the disc is shown, where the short 
dashed line indicates the one pressure scale height limit for the 
gas disc. The dash-dotted line represents the inclination of the 
binary orbit plane relative to the disc midplane. 
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collisions and measure relative velocities is identical to the 
one applied to the low inclination case (model 8). In FigfT2l 
(panel 1 and 2, solid lines) we present collisional velocities be- 
tween equally and differently sized planetesimals, respectively. 
Collisional velocities are even higher in this case (~ a few km/s) 
than for the low inclination case for all size combinations con- 
sidered. 

The orbital elements are displayed in FigfO] We observe 
that the Kozai effect causes an increase of eccentricity up to 
e = 0.4-0.6, such that the contribution to the relative velocity 
between planetesimals due to the term oc e ; Aw, becomes large, 
and this is the dominant cause of high velocity collisions be- 
tween equal sized bodies. The situation is different for planetes- 
imals of different sizes. The different gas drag strengths change 
the relative inclination between the planetesimals and the disc, 
resulting in a variation of the inclination perturbations caused 
by disc gravity, which on the one hand causes the value of 
that is approached during the first half of the Kozai cycle to be 
different, and on the other hand causes the Kozai mechanism 
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Fig. 12. Average collisional velocities in the high inclina- 
tion jp = 0.78(45°) case between equally (panel 1) and 
differently sized planetesimals (panel 2). Left panel: 10km- 
10km collisions (green-solid); lkm-lkm collisions (blue- 
solid); lOOm-lOOm collisions (red-solid). Right panel: 10km- 
lkm (green-solid); lOkm-lOOm (blue-solid); lkm-lOOm (red- 
solid). Threshold velocities for catastrophic disruption corre- 
sponding to the different size-combinations for the strong ag- 
gregates (dashed line) are also shown. Relative velocities have 
also been calculated using a larger orbit width Aa = 2 ■ 10~ 3 AU 
(dotted lines). Panels 3 and 4: collision count for orbital width 
Aa = 2 ■ 10~ 4 AU and Aa = 2 ■ 10~ 3 AU, respectively. The 
line colours and styles in panel 3 and 4 are as follows: green- 
solid (lOkm-lOkm); blue-solid (lkm-lkm); red-solid (100m- 
1 00m) ; green-dashed ( 1 0km- 1 km) ; blue-dashed ( 1 0km- 1 00m) ; 
red-dashed (1km- 100m). 
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Fig. 13. Average orbital parameters of the colliding planetes- 
imal pairs in the high inclination jp = 45° case. The line 
colours and styles are as follows: green-solid (lOkm-lOkm); 
blue-solid (lkm-lkm); red-solid (lOOm-lOOm); green-dashed 
(10km- lkm); blue-dashed (10km- 100m); red-dashed (1km- 
100m). 



to operate on different time scales. As a result planetesimals of 
different sizes decouple from the disc at different times. As can 
be seen from panels 4 and 5 of FigQj] the orbital planes and 
pericentres of different sized planetesimals become randomly 
distributed with A£2, ~ 1 and Aw, ~ 1, and the relative ve- 
locities are increased substantially due to the terms which are 
proportional to <?,Aw,, Aa,, and sin (ff,)Af2, in Eq. (l28l i, with the 
latter term being the dominant one. As found in the low incli- 
nation case, the relative velocities of different sized planetesi- 
mals are generated due to misalignment of their pericentres Aw, 
and their orbital planes AQ,. In contrast to the low inclination 
model, however, the eccentricities, <?,-, inclinations a,, and the 
differential nodal precession angle, A£2,, are much larger due to 
the Kozai effect, leading to very large collisional velocities. 

In FigQ~2](panel 3) we display the number of collisions de- 
tected for the different size combinations, and it may be ob- 
served that the number of collisions drops to ~ 100 for colli- 
sions between differently sized planetesimals. As the bodies be- 
gin to orbit in different planes their encounter proba bility is re- 
duced, as discussed in Sect. !5.2.T1 and reported by Kie & Zhou 
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d2009h . This implies that the averaged relative velocities (pan- 
els 1 and 2) are obtained from only ~ 10 2 — 10 3 data points, and 
the results are in danger of becoming statistically unreliable. In 
order to compare these collision velocities with more statisti- 
cally significant data, we also calculated the relative velocities 
while adopting a larger cross sectional area for the intersecting 
orbits corresponding to Ar = 2 • 10 3 AU. The collision prob- 
ability is increased by a factor of ~ 10, as observed in panel 4 
of FigQ~2] These relative velocities are plotted as dotted lines in 
panels 1 and 2 of Fig. [12] and are almost indistinguishable from 
the data obtained with the smaller orbit width, suggesting that 
the results are reliable. As was discussed for the low inclination 
model, the reduction in the collision frequency can have poten- 
tially important consequences for planetesimal accretion, since 
it may favour collisions between similar sized bodies which or- 
bit in similar planes. But, it is clear from FigJT2]fhat when the 
Kozai effect is active, large collision velocities are obtained for 
all size combinations. 

The disruption/erosion threshold velocity for the strong ag- 
gregates is shown using the dashed lines in panels 1 and 2. 
The collisional velocities are substantially larger than those re- 
quired for erosion or catastrophic disruption for all size combi- 
nations in this case, leading to the unsurprising conclusion that 
planetesimal accretion is not possible if the Kozai mechanism 
operates during the epoch of planet formation. For the relative 
velocities observed in FigfT2lof ~ 2km/s for equally sized plan- 
etesimals, we estimate that collisions will always lead to frag- 
mentation unless bodies of size ~ 10 3 km have already formed. 
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Fig. 14. Orbital elements for model 4 (left panels) and model 5 
(right panels). Only planetesimals with semi-major axes a < 13 
AU are included. In model 4 the Kozai effect is still operating, 
while in model 5 it is not. Panels 1 and 2 show the eccentricity. 
Panels 3 and 4 display the inclination, where the short dashed 
lines represent the inner and outer edge of the disc and the long 
dashed line indicates the threshold inclination of 0.68 (39.2°) 
above which the Kozai mechanism can operate. 



6. Suppressing the Kozai effect 

6.1. Increasing the disc mass (models 4 and 5) 

The Kozai effect generates relative velocities which are always 
too large for collisional growth of planetesimals, at least for 
sizes in the range (lOOm-lOkm). An interesting question is un- 
der what conditions does the Kozai mechanism cease to oper- 
ate. As discussed in Sect. 5.3 the Kozai mechanism relies on the 
apsidal precession induced by the binary companion halting, 
such that to, = 0, which then allows for a secular change of the 
eccentricity and inclination values. If apsidal precession can be 
induced by another source, however, such as the disk for exam- 
ple, then we would expect the Kozai mechanism t o be ineffec 



tive o nce this induced precession is fast enough ( Wu & Murray 
2003). 



Since the apsidal precession rate caused by the disc scales 
with disc mass (see Eqi30t. we have run simulations with an 
increased disc mass to see if the Kozai mechanism can be 
suppressed in this way. The results are shown in Fig JT4l for 
model 4 (left panels) and 5 (right panels) with disc masses of 
Md - 3.3M and 6.6M, respectively, where M is the nominal 
disc mass used in model 3. As can be seen from the figure, the 
Kozai mechanism still switches on for the Md = 3.3M case but 
does not operate in the Md = 6.6M model. Consequently the 
eccentricity grows to large values in model 4, as seen in panel 
1 of Fig [14] while it remains low in model 5 (panel 2). 



6.2. Increasing binary separation (models 6 and 7) 

Another way to increase the relative importance of the disc- 
induced apsidal precession is to decrease the gravitational in- 
fluence of the binary companion. In this section we present 
simulations in which the separation of the binary system was 
increased to examine when the Kozai effect is suppressed. We 
consider separations of D — 90 AU and D — 120 AU in models 
6 and 7, respectively. Because the Kozai time scale is oc D 3 , 
extremely long simulation times are required to prove the exis- 
tence (or non existence) of the Kozai effect. The required run 
times would be prohibitive for a full 3D hydrodynamic simula- 
tion, so we adopted the compromise of switching off the hydro- 
dynamic evolution of the disc in order to speed up the run times 
for these models. The disc in these simulations therefore re- 
mains axisymmetric and static in the precessing frame, but our 
adoption of a precessing reference frame causes it to precess 
around the binary angular momentum vector at the prescribed 
rate. 

We tested the accuracy of this approach by re-running 
model 3 using a static, precessing disc model, and found sim- 
ilar results when compared with the run in which the disc was 
allowed to evolve self-consistently. Differences between the 
static, precessing disc model and the full hydrodynamic sim- 
ulation arise mainly because the disc experiences a low am- 
plitude nodding motion (oscillation in its inclination) with a 
peri od equal to half the binary orbit period in the full simula- 
tion ( ILarwood et al .111 9961: iFraener & Nelsonkoiol) . where this 
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Fig. 15. Orbital elements for model 6 (left panels) and model 
7 (right panels). The Kozai effect is operating in model 6, but 
is not in model 7. The eccentricities are shown in panels 1 and 
2. The inclination is shown in panels 3 and 4. The long dashed 
lines in panels 3 and 4 represents the 39.2° inclination thresh- 
old, above which the Kozai effect can operate. 

is driven by the time-varying gravitational field of the binary. 
Although the effect is relatively small, it is likely to reduce 
the accuracy with which a rigidly precessing disc model can 
be used to determine planetesimal collision velocities in detail. 
But, it is a useful set-up for determining the disc mass for which 
the Kozai mechanism operates. 

The results are shown in Fig {15] for the D = 90 AU (left 
panels) and D = 120 AU model (right panels). We can observe 
that the Kozai effect is operating in model 6 while it is not in 
model 7. The time scale for the Kozai effect to increase the 
eccentricities is larger by a factor 3.3, as expected, for model 6 
compared to the reference model 3 due to the Kozai time scale 
increasing as ~ D 3 . If the Kozai effect was operating in model 
7, we would expect to see substantial eccentricity growth after 
500 orbits (a factor 8 longer than in our reference model). This 
is not seen, however, and we conclude that the Kozai effect 
cannot operate in model 7. 

6.3. A theoretical argument 

We have seen in the previous two sections that the Kozai mech- 
anism can be suppressed provided the disc mass or the binary 
separation are large enough. Now we want to understand the 
numerical results by means of a theoretical argument. As has 
already been discussed, the Kozai mechanism relies on halting 
the apsidal precession of the planetesimals induced by the bi- 
nary companion. If this can be prevented then there should be 
no secular net change of eccentricities or inclinations. The to- 
tal net rate of apsidal precession is given by the sum of the 
contributions by the binary companion (Eqj38l and the disc 



induced apsidal precession. Since we only consider cases for 
which 6i < aj, the latter is given by Eq. (l30l . Hence the to- 
tal apsidal precession rate experienced by a planetesimal on a 
circular orbit (e,- = 0) is: 



— - = —j (20a,-) i (2-5 sin 2 (w,) sin 2 (a,)) - \6j d \, 



dt 2D 3 



(44) 



where a and D are expressed in units of AU. The first term on 
the right hand side of Eq. (l44b accounts for the binary-induced 
precession, and (1>d is the disc-induced precession rate given 
by Eq. d30b . Both of these rates are expressed in time units of 
orbits at a = 20 AU (the time unit used throughout the paper). 
In order for the Kozai mechanism to operate we require that 
^ = 0. This can only be true if the prograde term due to the 
binary companion is larger than the retrograde term of the disk. 
For the Kozai effect to operate, we therefore require that: 



^ (20ad* > \0d\ 



(45) 



Using Eq. d30"i l. we can express this a condition on the total disc 
mass in terms of the planetesimal semi-major axis and binary 
orbit separation: 



3ttM 
M d < — (20a;) 5 



**«GSf)*«(sf) 



(46) 



If this inequality is fulfilled we expect the Kozai mechanism 
to operate. In Fig[16]we show how the upper limit on the disc 
mass varies as a function of binary separation in order for the 
Kozai mechanism to just switch on/off for a planetesimal or- 
biting at a representative value of the semi-major axis a,- = 11 
AU, which is approximately half of the disc tidal truncation ra- 
dius. The red area corresponds to the regime where the inequal- 
ity is fulfilled and the Kozai-mechanism should operate. The 
cyan area marks the parameter space where the Kozai mecha- 
nism should not operate (large disc mass/large binary separa- 
tion). The symbols represent the simulation results, where open 
squares indicate that the Kozai mechanism was found to oper- 
ate (models 3, 4 and 6), while crosses represent the models in 
which the Kozai effect was ineffective (models 5 and 7). We 
see that the numerical results agree with our predictions. 

We also plot the boundary lines that correspond to semi- 
major axes of a, = 6 AU (dashed line) and a, = 15 AU (dot- 
ted line). Between a, = 6 AU and a, ■ - 11 AU the precession 
frequency due to the disc is nearly independent of semi-major 
axis, as can be seen from the solid line in panel 7 of FigQ] Yet 
the binary precession still increases as oc a 2 . Hence the param- 
eter space for which the Kozai mechanism operates becomes 
larger as the semi-major axis increases. Between a, ■ — 1 1 AU 
and a, = 15 AU the disc precession rate increases in magni- 
tude approximately as oc a. . Thus the boundary separating 
the Kozai-active versus Kozai-inactive regions does not change 
very much beyond a, ■ - 11 AU. 

7. Discussion and conclusions 

In this paper we have investigated the dynamics of planetesi- 
mals that are embedded in a gaseous disc that is perturbed by 
a binary companion on a circular, inclined orbit. In contrast to 



22 



M.M. Fragner, R.P. Nelson & W. Kley: On collisional growth of planetesimals in misaligned binary systems 



Disk mass versus binary separation 



,_, 6 




Fig. 16. Plot of the parameter regime explored in the simu- 
lations. The diagram shows disc mass against binary separa- 
tion, with the red area denoting the region where the Kozai 
mechanism should operate, and the cyan area showing the re- 
gion where it should not for a planetesimal located at 1 1 AU. 
The symbols denote the outcome of the simulations, where a 
cross implies that the Kozai effect was switched off, and an 
open square indicates that the Kozai mechanism operates. The 
dashed and dotted line represents the same boundary for a body 
at a,- = 6 AU and a-, — 15 AU respectively. 



previous work the planetesimals are allowed to interact with 
the gaseous disc via gas drag and disc gravity. The disc evolves 
hydrodynamically because of gravitational perturbations due to 
the binary companion, and the disc model we consider under- 
goes solid body precession around the orbital angular momen- 
tum vector of the binary system, as expected when the warp 
propagation time via bending waves is shorter than the differ- 
ential precession time. Discs around young stars are expected 
to show similar behaviour since it is estimated that a < H/R 
in these systems, where a is the usual viscosity parameter. The 
key results of our study can be summarized as follows. 

- In the absence of disc gravity, planetesimals undergo strong 
differential nodal prece ssion, such that they eventually or- 



bit in different pl anes (IMarzari. Thebault & Scholll 2009 



Xie & Zhou 2009). They also precess relative to the disc, 



so their orbits become highly inclined relative to it. We find 
that disc gravity acts to prevent this differential precession, 
and forces the planetesimals to precess with the disc on av- 
erage. Viewed locally, the forcing of nodal precession by 
the binary companion causes the planetesimals to oscillate 
about the midplane of the disc, with an amplitude that de- 
pends on the size of the planetesimals because of the influ- 
ence of gas drag in damping this relative motion. This key 
result suggests that the influence of the gravitational field 



of the precessing disc should always be included in future 
studies of planet formation in inclined binary systems. 

- Previous studies that focused on the influence of a coplanar, 
eccentric binary companion have shown that the eccentric 
binary generates a forced eccentricity in the planetesimal 
swarm. Gas drag acts to align the pericentres of planetes- 
imals which are of the same size, such that collisions be- 
tween same-sized bodies are dominated by keplerian shear. 
Collisions between different sized bodies occur with higher 
veloc i ty due to the misaligned orbits (IMarzari & Scholl 
2000; iThebault. Marzari & Scholll 120061) . We find that for 
a circular binary, where the eccentricity of planetesimals is 
largely generated by high frequency terms in the disturbing 
function, strong alignment of pericentres is not observed. 

- Previous work which has examined planetesimal dynamics 
in modestly inclined and eccentric binary systems suggests 
that gas drag causes size-dependent phasing of both peri- 
centres and the lin es of nodes of perturbed planetesimals 
JXie & Zhoul l2009). This has the effect of favouring colli- 
sions between same sized objects, for which the collision 
velocity is dominated by the keplerian shear, and as such it 
has been suggested that this may provide an effective chan- 
nel for planetesimal growth in inclined binary systems. Our 
results obtained for a binary with inclination yp = 25° are 
in basic agreement with this, as we find that same-sized 
planetesimals do indeed occupy very similar orbital planes, 
whereas the orbits of differently sized bodies are mutually 
inclined. This arises in our case because of the different 
amplitudes of oscillation about the disc midplane observed 
for planetesimals of different size, with smaller bodies re- 
maining closer to the midplane. We thus find that the col- 
lision velocities between differently sized bodies are sig- 
nificantly larger than between same-sized bodies (typically 
Vcoii — 200 ms _1 for different sized objects and v co ii — 50 

- 70 ms _1 for same-sized objects). We note that the lack of 
pericentre alignment observed in our simulations leads to 
the substantial collision velocities between same size bod- 
ies. Collisions with the velocities described above are likely 
to be erosive or disruptive, depending on the mass ratio of 
the colliding bodies. 

- For highly inclined systems with yp = 45°, we find that 
the Kozai mechanism can operate, causing large changes in 
the eccentricities and inclinations of the planetesimals, and 
leading to collisional velocities that are much too large to 
allow for planetesimal accretion. The long term influence 
of the disc in cases where the Kozai mechanism operates 
is to modify the period associated with the Kozai cycle by 
periodically forcing the planetesimal inclinations to fall be- 
low, or rise above, the critical value for the Kozai effect to 
operate, an = 39.2°. We find that planetesimals of size 10m 

- 10km all experience the Kozai effect. 

- Increasing the disc mass or binary separation can suppress 
the Kozai effect, as disc gravity can induce apsidal preces- 
sion which is fast enough to render the Kozai mechanism 
ineffective. For a binary separation of D — 60 AU, and disc 
truncation radius of 20 AU, we find that we need to in- 
crease the disc mass by a factor of ~ 6 above the minimum 
mass solar nebula (MMSN) value in order to switch off the 
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Kozai effect at a representative orbital radius of a, — 11 AU. 
Increasing the binary separation to D — 120 AU, and main- 
taining the mass of our disc at its nominal value (equivalent 
to the MMSN) also rendered the Kozai mechanism inoper- 
ative. 

In the light of the above findings, we can conclude that 
highly inclined, distant binary companions probably will not 
induce the Kozai effect during planet formation while the disc 
is present, although it may do so once the disc has dissipated, 
such that the final planetary system is strongly perturbed. But 
we note that the peak in the distri bution of binary orbit sep 



arations occurs a t D ^ 30 AU ([D uquennov & Mayor] [1991 



Ghe z et alJl!993t iLeinert et al.l ll993). which is small enough 



for the Kozai mechanism to operate during the planet formation 
epoch when the gas disc is present. There are number of extra- 
solar planetary systems, however, which are observed to be in 
relatively close binary systems, y Cephei being a notable exam- 
ple (see Kley & Nelson 2008, and references therein). Here the 
binary separation is only a few tens of AU, such that the Kozai 
mechanism could possibly have operated when planetesimals 
were accreting. The fact that a planet is observed in this system 
with a modest orbital eccentricity is probably an indication that 
the binary inclination is too small to allow the Kozai m echa- 
nism t o operate. This is consistent with the findings of Hale 
(1 1994b that binary orbits with D < 40 AU tend to be reason- 
ably well aligned with the stellar spin axes. 

A key issue that needs to be addressed in the future is the 
evolution of a planetesimal swarm in a binary system where the 
orbit is both inclined and eccentric. Our results suggest that a 
circular binary system does not allow for the pericentres of per- 
turbed planetesimals to be well aligned, since a circular orbit is 
unable to impose a preferred direction on the system, and this is 
in clear contrast to studies which have considered an eccentric 
companion. It should be noted, however, that the inclusion of 
the disc gravity will be of crucial importance in such a study, 
since it will significantly perturb the orbits of the planetesi- 
mals. This is particularly the case when the planetesimal orbits 
are inclined relative to the disc midplane, and the presence of 
an eccentric binary c ompanion may also cause the disc itself 



to become eccentric dKlev & Nelson! 120081: iKlev etalj [2008) 



(where the disc eccentricity may be dependent on the disc self- 
gravity (IMarzari et al. 2009)), further complicating matters. We 
further note that the inner edge of our computational domain 
was located at R - 2 AU, so our present study examines planet 
formation in the region normally associated with giant planet 
formation. Studies of this type which examine terrestrial planet 
formation closer to the central star will require deployment of 
substantially more powerful computational resources than were 
used here, in order to simulate the larger range of time scales in 
the problem. Such a study should also include a more realistic 
planetesimal size distribution so that the possibilty of accretion 
occuring via collisions between similarly sized bodies can be 
explored. 

In this paper we have only considered planet formation 
via collisions between planetesimals, whereas dust accumu- 
lation onto individual planetesimals could still lead to plane- 
tary growth as long as collisions between planetesimals can 



be avoided. This issue has been r ecently been ad dressed 
bv lPaardekooper & Leinhardtl (I2010l) : IXie et al.l (l2010l) . where 
they show that km-sized planetesimals may grow in size by 
two orders of magnitude due to accretion of collisional debris 
if the efficiency of planetesimal formation stays low, i.e. only a 
few planetesimals form. This study involved two-dimensional 
simulations, and it will be interesting to examine how the re- 
sults change for a misaligned system where planetesimals may 
orbit out of the disc midplane where the majority of the dust 
and collisional debris will reside. As we have seen in Sect. 5.2, 
10 km sized planetesimals tend to have orbits which are in- 
clined relative to the disc by more that the disc scale height, 
once their free particle rates are substantially different from the 
disc precession rate. Planetesimal formation via dust accretion 
is therefore only expected in radial regions of the disc where 
the precession rate of the planetesimals and the disc are well- 
matched. Closer to the inner or outer edge of the disc, efficient 
dust accretion may not be possible due to the large relative in- 
clinations of the planetesimals. This, and other issues discussed 
in this paper, will be the subject of future publications. 
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